Kuiper Belt-Like Hot and Cold Populations of Planetesimal Inclinations in the $\beta$ Pictoris Belt Revealed by ALMA

TL;DR

This study uses ALMA observations to reveal that the $eta$ Pictoris planetesimal belt has a bimodal inclination distribution similar to the Kuiper belt, indicating complex dynamical history involving hot and cold populations.

Contribution

First direct imaging evidence of a bimodal inclination distribution in an extrasolar planetesimal belt, suggesting multiple dynamical populations and challenging simple formation models.

Findings

Inclination distribution is better described by two Gaussian populations.

The hot and cold populations have RMS inclinations of approximately 9 and 1 degree.

The belt's midplane is misaligned with $eta$ Pic b's orbital plane.

Abstract

The inclination distribution of the Kuiper belt provides unique constraints on its origin and dynamical evolution, motivating vertically resolved observations of extrasolar planetesimal belts. We present ALMA observations of millimeter emission in the near edge-on planetesimal belt around $\beta$ Pictoris, finding that the vertical distribution is significantly better described by the sum of two Gaussians compared to a single Gaussian. This indicates that, as for the Kuiper belt, the inclination distribution of $\beta$ Pic's belt is better described by the sum of dynamically hot and cold populations rather than a single component. The hot and cold populations have RMS inclinations of 8.9$^{+0.7}_{-0.5}$ and 1.1$^{+0.5}_{-0.5}$ degrees. We also report that an axisymmetric belt model provides a good fit to new and archival ALMA visibilities, and confirm that the midplane is misaligned with respect to $\beta$ Pic b's orbital plane. However, we find no significant evidence for either the inner disk tilt observed in scattered light and CO emission or the South-West/North-East (SW/NE) asymmetry previously reported for millimeter emission. Finally, we consider the origin of the belt's inclination distribution. Secular perturbations from $\beta$ Pic b are unlikely to provide sufficient dynamical heating to explain the hot population throughout the belt's radial extent, and viscous stirring from large bodies within the belt alone cannot reproduce the two populations observed. This argues for an alternative or additional scenario, such as planetesimals being born with high inclinations, or the presence of a `$\beta$ Pic c' planet, potentially migrating outwards near the belt's inner edge.

Figures

Tables

Table 1. Table 1: Best-fit model parameters

Parameter	Unit	Model 1	Model 2	Model 3	Model 4	Model 5
$r_c$	au	${104.9}_{-1.1}^{+1.1}$	${107.5}_{-1.0}^{+1.0}$	${106.3}_{-1.3}^{+1.4}$	${107.8}_{-1.0}^{+1.0}$	${105.2}_{-1.4}^{+1.4}$
$\sigma$	au	${38.9}_{-1.1}^{+1.3}$	${36.7}_{-1.3}^{+1.0}$	${36.6}_{-1.2}^{+1.1}$	${36.4}_{-1.3}^{+1.1}$	${36.3}_{-1.2}^{+1.1}$
$h$		${0.070}_{-0.004}^{+0.004}$	${0.014}_{-0.006}^{+0.006}$	${}_{}{}^{\ast }0.020_{-0.008}^{+0.010}$	${0.054}_{-0.004}^{+0.004}$	${}_{}{}^{\ast }0.072_{-0.012}^{+0.014}$
$I$	∘	${86.6}_{-0.3}^{+0.4}$	${88.6}_{-0.3}^{+0.2}$	${88.8}_{-0.3}^{+0.4}$	${89.3}_{-0.5}^{+0.4}$	${89.5}_{-0.4}^{+0.3}$
$PA$	∘	${30.0}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$	${29.8}_{-0.1}^{+0.1}$	${29.7}_{-0.1}^{+0.1}$
$h_2$		-	${0.110}_{-0.006}^{+0.008}$	${}_{}{}^{\ast }0.14_{-0.02}^{+0.02}$	-	-
$a$		-	${0.20}_{-0.04}^{+0.05}$	${0.21}_{-0.04}^{+0.05}$	-	-
$\beta$		-	-	${0.7}_{-0.2}^{+0.2}$	-	${0.4}_{-0.2}^{+0.2}$
$p$		-	-	-	${0.9}_{-0.1}^{+0.1}$	${0.85}_{-0.07}^{+0.07}$
Free Parameters		14	16	17	15	16
${\chi }^2$		1456364.9	1456297.6	1456297.4	1456312.2	1456297.3
$\Delta {\chi }^2$		-	-67.3aaDifference with respect to Model 1	-0.2bbDifference with respect to Model 2	+14.6bbDifference with respect to Model 2	-0.3bbDifference with respect to Model 2
$\Delta$AIC		-	-63.3aaDifference with respect to Model 1	+1.8bbDifference with respect to Model 2	+12.6bbDifference with respect to Model 2	-0.3bbDifference with respect to Model 2
$\Delta$BIC		-	-38.9aaDifference with respect to Model 1	+14.0bbDifference with respect to Model 2	+0.33bbDifference with respect to Model 2	-0.3bbDifference with respect to Model 2

Equations

$$\rho ={\rho }_0{e}^{-\frac{(r-r_{\mathrm{c}}{)}^2}{2{\sigma }^2}}\frac{e^{-\frac{z^2}{2(hr{)}^2}}}{\sqrt{2\pi }hr}$$

$$\rho ={\rho }_0{e}^{-\frac{(r-r_{\mathrm{c}}{)}^2}{2{\sigma }^2}}\left[\frac{a}{\sqrt{2\pi }hr}{e}^{-\frac{z^2}{2(hr{)}^2}}+\frac{1-a}{\sqrt{2\pi }{h}_{2}r}{e}^{-\frac{z^2}{2(h_{2}r{)}^2}}\right].$$

$$\rho ={\rho }_0{e}^{-\frac{(r-r_{\mathrm{c}}{)}^2}{2{\sigma }^2}}C{e}^{-{\left(\frac{|z|}{\sqrt{2}hr}\right)}^p},$$

$$L(\theta )={\int }_{\theta }^{\pi /2}f(i)\frac{\cos \theta }{\sqrt{{\sin }^{2}i-{\sin }^2\theta }}di.$$

$$L(\theta )\sim \frac{\cos \theta }{{\sigma }_i^2}{\int }_{\theta }^{\pi /2}\frac{i{e}^{-\frac{i^2}{2{\sigma }_i^2}}}{\sqrt{i^2-{\theta }^2}}di,$$

$$L(\theta )=\sqrt{\frac{\pi }{2{\sigma }_i^2}}\cos \theta \mathrm{erf}\left(\sqrt{\frac{\frac{{\pi }^2}{4}-{\theta }^2}{2{\sigma }_i^2}}\right)e^{-\frac{{\theta }^2}{2{\sigma }_i^2}}.$$

$$\rho (z)\propto e^{-\frac{z^2}{2{\sigma }_z^2}},$$

$$\bar{i}=\sqrt{\frac{\pi }{2}}h\mathrm{and}{\mathrm{RMS}}_i=\sqrt{2}h.$$

$$\cos (\mathrm{d}I)=\cos (I_{\mathrm{b}})\cos (I_{\mathrm{p}})+\sin (I_{\mathrm{b}})\sin (I_{\mathrm{p}})\cos ({\Omega }_{\mathrm{p}}-{\Omega }_{\mathrm{b}}).$$

$$v_{\mathrm{rel}}=\sqrt{6⟨i^2⟩}{v}_{\mathrm{Kep}}=103.9M_{\star }^{0.5}h{r}^{-0.5},$$

$$\sqrt{⟨i^2⟩}={\left(2C_{\mathrm{I}}\frac{M}{M_{\star }^2}{r}^2\Omega \Sigma t\right)}^{\frac{1}{4}},$$

$$M\Sigma =4.43\times 10^3{M}_{\star }^{\frac{3}{2}}{r}^{-\frac{1}{2}}\frac{{\sqrt{⟨i^2⟩}}^4}{t}{\mathrm{M}}_{\oplus }^2{\mathrm{au}}^{-2},$$

$$R_{\mathrm{col}}=\sigma n{v}_{\mathrm{rel}}\left(1+\frac{v_{\mathrm{esc}}^2}{v_{\mathrm{rel}}^2}\right),$$

$$\Sigma >\frac{M}{2\pi r\Delta r_{\mathrm{stir}}},$$

Full text

submillimetre: planetary systems – planetary systems – circumstellar matter – Kuiper belt: general – protoplanetary disks – stars: individual ($\beta$ Pictoris (catalog )).

1 Introduction

The spatial distribution of dust in planetesimal belts around main-sequence stars provides significant constraints on the formation and dynamical evolution of planetary systems, both individually and as a population (e.g. Marino et al., 2018; Matrà et al., 2018a). Due to its vicinity and youth, the A6-type main sequence star $\beta$ Pictoris has long been considered a unique laboratory to understand the outcome of the planet formation process. Over more than three decades since the discovery of an infrared excess above the stellar emission (Aumann, 1985), indicative of the presence of a dusty disk, numerous studies have been attempting to unravel the complexity of this planetary system.

At only 19.44 pc (van Leeuwen, 2007), the disk’s brightness and on-sky extent have made this system extremely suitable for resolved imaging. It was as early as 1984 when the belt’s near edge-on geometry was discovered through optical coronagraphic observations (Smith & Terrile, 1984), which showed starlight scattered off micron sized particles extending out to beyond 1000 AU. Soon after, improved scattered light imaging of this system from the ground at optical and near-infrared (near-IR) wavelengths unveiled a complex disk structure with several asymmetries (Kalas & Jewitt, 1995).

The near edge-on configuration makes $\beta$ Pic’s belt extremely well suited for studies of its vertical structure. Amongst the asymmetries, a warp in the spine of the on-sky brightness distribution within $\sim 80$ au from the star was first unveiled through HST imaging and subsequently confirmed in all follow-up optical/near-IR scattered light observations (Mouillet et al., 1997; Heap et al., 2000; Golimowski et al., 2006; Boccaletti et al., 2009; Lagrange et al., 2012; Milli et al., 2014; Apai et al., 2015; Millar-Blanchaer et al., 2015). The presence of this warp was soon attributed to gravitational interaction with an unseen companion in the inner regions of the system (Mouillet et al., 1997; Augereau et al., 2001). It wasn’t until 2009 that a giant planet, $\beta$ Pictoris b, was discovered through direct imaging observations (Lagrange et al., 2009, 2010) at a projected distance of $\sim$9 AU from the star (Lagrange et al., 2018a), misaligned to the outer belt (Lagrange et al., 2012) in a way that is roughly consistent with the general expectation from the observed scattered light warp (e.g. Dawson et al., 2011; Nesvold & Kuchner, 2015). Interestingly, the position angle of the CO and C I gas disks is similar to that of the inner disk seen in scattered light, suggesting a link between the two (Matrà et al., 2017; Cataldi et al., 2018), but no constraints have been reported so far on the presence of a warp in dust emission at longer wavelengths, due to a lack of angular resolution necessary to resolve this feature.

So far, the observed vertical distribution and width of the belt has received less attention. HST observations found that a sum of two Gaussians or Lorentzians best describes the observed vertical distribution of scattered light emission (e.g. Golimowski et al., 2006; Lagrange et al., 2012; Apai et al., 2015), with the scale height of the broad Gaussian component reaching values as high as $\sim$12 au at a distance of $\sim$80 au from the star. This corresponds to an aspect ratio of $\sim 0.15$ at 80 au, similar to the aspect ratio inferred from modelling of polarimetric near-IR observations ($\sim$0.137, Millar-Blanchaer et al., 2015).

The vertical distribution of planetesimal belts is a crucial observable as it reflects the inclination distribution of the constituent particles. In our own Kuiper belt, the inclination distribution is broad and bimodal (e.g. Brown, 2001; Kavelaars et al., 2009; Petit et al., 2011), and sets tight constraints on its dynamical history (Morbidelli et al., 2008), particularly on the outward migration of Neptune in the early stages of the Solar System’s evolution (e.g. Nesvorný, 2015). Therefore, accessing the inclination distribution of particles tracing planetesimals in extrasolar planetesimal belts is an attractive avenue to probe their dynamics and the effects of unseen planets.

With the assumption of an externally unperturbed and initially cold belt, the inclination distribution (particularly its width) could for example yield a direct link to the relative velocities of planetesimals. This allows us to probe the level of dynamical heating induced by unseen, large gravitational perturbers within the belt (e.g. Quillen et al., 2007). To date, only observations in optical/near-IR scattered light have had sufficient resolution to vertically resolve edge-on belts such as around $\beta$ Pic. Unfortunately, however, the radiation pressure felt by the small grains probed by these observations can significantly excite their eccentricities leading to a consequential increase in their inclinations (Thébault, 2009), which has prevented drawing a robust link between the observed vertical structure and the presence of large, unseen bodies required to create it.

In this paper, we tackle this issue using new ALMA interferometric observations that vertically resolve the $\beta$ Pictoris belt at 1.33 mm. The observations are described in §2, and the main results derived from image analysis in §3. In §4, we confirm these results by fitting the interferometric visibilities using parametric models to describe the belt’s dust density distribution. In §5, we link the belt’s observed vertical distribution to the inclination distribution of mm grains, and discuss its origin in the context of gravitational perturbations by large bodies within or exterior to the belt. Finally, in §6, we draw our conclusions and summarize our findings.

2 Observations

ALMA observations of the $\beta$ Pictoris disk were obtained within its Cycle 2 (project code 2012.1.00142.S) using Band 6 receivers. Data was taken with the 12-m array in extended (one pointing centered on the star) and compact configuration (two pointings centered $\pm 5\arcsec$ from the star along the belt midplane), as well as with the Atacama Compact Array (ACA, one pointing centered on the star). Baselines covered a range between 9 and 1574 m, enabling us to be sensitive to structure between about 0$\farcs$3 and 27$\arcsec$. A more complete description of the observations, including calibration, reduction and imaging of CO and other lines can be found in Matrà et al. (2017, 2018b). For continuum imaging, we flagged spectral channels around the frequency of the detected CO J=2-1 line (230.538 GHz), and combined all the datasets from the different ALMA configurations as well as the ACA. In doing so, we applied a constant rescaling factor to the visibility weights of each dataset, where these factors (different for each dataset) were determined from the scatter of the observed visibilities in §4.1. Then, we imaged the combined dataset using the CLEAN algorithm (Högbom, 1974) in multi-frequency synthesis mode. The dataset covers $\sim$6.8 GHz of total bandwidth at an effective wavelength (frequency) of 1.33 mm (224.58 GHz). We used natural weighting of the visibilities to ensure maximum sensitivity while achieving a synthesized beam of size $0\farcs 32\times 0\farcs 27$ and PA of -82$\fdg$01; this implies a resolution of $6.1\times 5.3$ au at the distance of the star (19.44 pc).

3 Results

3.1 New 1.33 mm dust continuum dataset

Fig. 1 shows the continuum emission map of the $\beta$ Pictoris system obtained from the combined, naturally weighted, visibility dataset. A peak signal-to-noise ratio (S/N) of 16 is reached at the location of the central star, for the RMS noise level of 12 $\mu$Jy/beam measured in an emission-free region of the continuum map. Emission from the belt is also clearly detected with a peak S/N of $\sim$14, as measured to the SW of the star along the midplane. The spatially integrated 1.33 mm flux measured on the primary-beam-corrected map within a $16\arcsec\times 4\arcsec$ box centered on the star is 20$\pm$2 mJy (where the uncertainty is dominated by the 10% contribution from flux calibration, Fomalont et al., 2014).

The stellar contribution as measured from the peak flux at its location is 190 $\mu$Jy, which should strictly be considered an upper limit due to extra disk emission lying along the line of sight to the star. In order to separate the stellar flux from the disk emission, we extract the subset of the dataset containing visibilities at $u-v$ distances $>300k\lambda$. This $u-v$ cutoff was chosen through inspection of CLEAN images, as an optimal trade-off to both ensure that all large-scale emission from the disk is filtered out, while retaining as much data as possible for maximum sensitivity. The obtained CLEAN image has an RMS noise level of 18 $\mu$Jy/beam for a beam size of $0\farcs 22\times 0\farcs 19$, and a compact point source is detected at the stellar position. A simple 2D Gaussian fit within CASA yields a stellar flux of 80$\pm$20 $\mu$Jy at 1.33 mm. This detection is consistent with the expected stellar flux of $\sim$86 $\mu$Jy obtained by extrapolating the Rayleigh-Jeans tail of the stellar blackbody-like emission from near-IR wavelengths. We therefore assume a stellar flux of 86 $\mu$Jy at 1.33 mm throughout the rest of this paper.

3.2 Belt radial structure

In Fig. 2, we show the profile of the continuum emission vertically integrated within $\pm 20$ au of the disk midplane, defined here as the position angle of the main disk observed in scattered light (29.3*∘*, Lagrange et al., 2012). Contrary to what is seen for CO emission in the same dataset (green line), the millimetre-sized dust traced by the 1.33 mm continuum observations (blue line) shows no significant peak brightness asymmetry between the SW and NE sides of the disk. This is in contrast with the SW/NE brightness asymmetry, between 30 and 80 au on either side the star, reported at 885 $\mu$m by Dent et al. (2014). The SW/NE ratio measured within the same radial distances as Dent et al. (2014), vertically integrated between $\pm 20$ au, is 1.03$\pm$0.04 in the new 1.33 mm dataset, consistent with no asymmetry. The same SW/NE ratio is 1.05$\pm$0.03 in the 885 $\mu$m archival dataset, indicating again a negligible (5$\pm$3)% asymmetry in this radial region. As this measurement can be affected by imaging in the presence of incomplete $u-v$ sampling and subsequent deconvolution via the CLEAN algorithm, we postpone further discussion of this issue in the context of visibility models to §4.6.

3.3 Belt vertical structure

Fig. 3 (top) shows vertical profiles of emission perpendicular to the disk midplane as a function of distance to the star, which were obtained by integrating emission every 15 au along the radial (midplane) direction. As shown by comparing the profiles (blue shaded region) to the resolution element of our observation (black dotted line), the dust disk is clearly resolved in the vertical direction. Following the procedure described in detail in §3.2 and Appendix B of Matrà et al. (2017), we fit Gaussians (purple dashed lines) to the vertical profiles at each midplane location, and derive vertical best-fit Gaussian centroids (forming the disk spine, Fig. 3 middle panel) and widths (measured as FWHM, Fig. 3 bottom panel). We report two main findings:

In contrast with the CO J=2-1 line emission, we find that dust emission is not well represented by a single vertical Gaussian, being more centrally peaked and having broader wings extending out to $\pm$40 au from the sky-projected midplane. This is evident in Fig. 4, which shows the residual vertical profile obtained by integrating emission across the entire midplane after vertically centering it so that the centroid of the emission aligns with the midplane. The best-fit Gaussian profile (red, with a standard deviation of $8.2\pm 0.2$ au) leaves significant residuals (Fig. 4, bottom), whereas a double-Gaussian profile (green, with standard deviations of $5.1\pm 0.3$ and $15.7\pm 1.4$ au) considerably improves the fit, leaving no significant residual emission. We find that relative flux calibration systematics of $\pm 10\%$ between the datasets with compact and extended baselines (tracing vertically broader and narrower emission, respectively) cannot account for this discrepancy, which must therefore be physical in origin. This shows that the double-Gaussian or Lorentzian profiles needed to fit scattered light observations (Golimowski et al., 2006; Lagrange et al., 2012; Apai et al., 2015) reflect the vertical structure of the millimeter emission that traces the dust-producing planetesimals. 2. 2.

The belt spine (Fig. 3, middle) presents little vertical displacement along the midplane and is largely consistent with the PA of the main disk observed in scattered light. This is significantly different from the CO J=2-1 emission, whose spine presents an extra tilt angle dPA of $\sim$4*∘* similar to that of the warp interior to $\sim$80 au seen in scattered light (Apai et al., 2015; Matrà et al., 2017). Therefore, we conclude that an inner disk tilt as large as observed in scattered light and CO observations is not detected in the structure of millimetre grains. Assuming the mm grains are tracing the planetesimal distribution, this shows that planetesimals are not perturbed into a warp - or at least not into one as pronounced as that observed in the distribution of the smallest grains observed in scattered light. At the same time, the tilt present in both the scattered light and the CO distribution, with the CO clump clearly offset by 4-6 au above the main disk midplane, suggests a connection between the two.

4 Modelling

In this section, we aim to confirm our results inferred from images through radiative transfer modelling of the interferometric visibility datasets, enabling us to constrain the 3D distribution of mm dust in the belt while avoiding potential artifacts introduced by the imaging and CLEAN deconvolution processes.

4.1 Method

We model the $\beta$ Pic belt using mm dust density distributions described in the following subsections, and the star as a point source with a flux density of 86 $\mu$Jy located in the geometrical center of the circular belt. For each model realization, we solve the radiative transfer at 1.33mm using RADMC-3D111http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/, and assuming the dust radial temperature profile expected from a blackbody around a star of 8.7 L*⊙* (Kennedy & Wyatt, 2014). Since we are interested in the structure rather than the overall amount of material in the belt, we fix the dust opacity and fit for the belt mass as a normalization factor, noting that this does not affect our conclusions on the belt structure. In addition to the parameters that govern the belt dust density distribution and its overall mass, we fit for belt geometry parameters (inclination $I$ and position angle PA).

Having created the model belt image, we produce visibility data separately for each of the ALMA arrays (12m and ACA), for each configuration of the 12m array and for each of the two mosaic pointings used for the compact configuration observations. We do so by first shifting the model belt image by RA and Dec offsets that we leave as free parameters in the fit, and multiplying it by the primary beam response of the relevant array, at its appropriate sky location for each observation. We then Fourier transform the image to produce a grid of model complex visibilities, which we evaluate at the $u-v$ locations sampled by the observation in question as described e.g. in Marino et al. (2016) and Tazzari et al. (2018). Finally, we add a point source representing the star to the model visibilities, located in the geometric center of the belt, and with its flux appropriately attenuated given its location with respect to the primary beam center.

The model visibilities are then fitted to the observed visibilities through a Monte Carlo Markov Chain (MCMC) method implemented through the emcee package (Foreman-Mackey et al., 2013). Eight nuisance parameters for the spatial shifts (in RA and Dec, so two for each dataset), two parameters to determine the viewing geometry ($I$ and PA), one normalization parameter setting the mass of the belt for a fixed opacity, and n parameters governing the dust distribution make up the 11+n - dimensional parameter space that we explore in our fit. Through the emcee package, we sample the 11+n dimensional posterior probability distribution of the explored parameters using the affine-invariant MCMC ensemble sampler of Goodman & Weare (2010). The posterior probability distribution is computed assuming uniform priors on all of the probed parameters, and a likelihood function proportional to $e^{-\chi^{2}/2}$, where $\chi^{2}$ is the usual definition of the chi-square function. The MCMC was run for each tested model with a number of walkers equal to 10x the number of free parameters, and over 2500 steps, ensuring each of the chains had reached convergence.

The uncertainty $\sigma$ on each visibility datapoint is defined using the visibility weights $w$ delivered by the ALMA observatory, where $\sigma_{u_{i}-v_{i}}=1/\sqrt{w_{u_{i}-v_{i}}}$. While assuming the delivered weights and hence the uncertainties are correct for different $u-v$ points relative to each other, we left a constant rescaling factor, equal for all $u-v$ points within a given dataset, as a free parameter in our model fit of §4.2. This is justified by previous ALMA modelling studies also needing to rescale the weights to appropriately represent the true scatter in the visibilities (e.g. Marino et al., 2018). We find this rescaling factor to be different but very well constrained to $0.4508\pm 0.0006$, $0.927\pm 0.014$, $0.690\pm 0.003$, and $0.698\pm 0.003$ for the 12m extended, ACA and for the two pointings of the 12m compact configuration, respectively. We therefore fix the rescaling factors to these values before producing the image in Fig. 1, and for all subsequent model fitting runs.

4.2 Radially and vertically Gaussian, axisymmetric model with radially constant aspect ratio

We begin by assuming the simplest dust density distribution model with the minimum number of free parameters, namely a radially and vertically Gaussian, axisymmetric belt with radially constant aspect ratio. The distribution of the dust density $\rho$ in the belt reads:

[TABLE]

where $r$ and $z$ are cylindrical coordinates, $r_{\rm c}$ and $\sigma$ are the center and standard deviation of the radially Gaussian belt, $h$ is the aspect ratio (here independent of radius), and $\rho_{0}$ is a normalization factor proportional to the total dust mass in the belt, which we fit for.

The best-fit parameters of this model, defined as the $50^{+34}_{-34}$th percentiles of their marginalized posterior probability distributions, are listed in Table 1 (Model 1). Within the framework of this model, and in qualitative agreement with Fig. 1, we find the disk to be near edge-on, with an inclination from face-on of $86.6^{+0.4}_{-0.3}$ degrees, and with its Gaussian radial distribution being centred at $\sim$105 AU with a width of $\sim$92 AU.

Fig. 5, top left shows the residuals obtained after subtracting the best-fit model visibilities from the data, and imaging the residual visibility dataset. The bottom left panel of Fig. 5 shows the residual radial profile constructed as in Fig. 2; the lack of significant residuals and SW/NE asymmetries reaffirms our conclusion that no clear departure from axisymmetry is present, within our measurement uncertainties, in the ALMA 1.33 mm dataset.

On the other hand, the right panel shows the vertical profile obtained after radially integrating the belt’s residual emission as done in Fig. 4. This displays significant residuals very similar to those obtained by simply subtracting a single Gaussian from the vertical flux distribution in the image plane. This confirms that the vertical distribution of mm dust in the planetesimal belt around $\beta$ Pictoris cannot be reproduced by a single Gaussian.

4.3 Double-Gaussian vertical structure model, with radially constant aspect ratio

Motivated by the appearance of the vertical emission distribution in the image plane, we construct a belt model with a vertical distribution characterised by a narrower and broader Gaussian. The Gaussians are both centered in the belt’s midplane ($z=0$) and are characterized, respectively, by aspect ratios $h$ and $h_{2}$ leading to standard deviations (scale heights) $hr$ and $h_{2}r$ at any given radius $r$. They comprise a fraction $a$ and $1-a$ of the total belt mass. The functional form of the dust density distribution thus reads

[TABLE]

This model is significantly better at describing the vertical distribution of the belt, as seen when comparing its residuals (Fig. 6) to the single vertical Gaussian model (Fig. 5). Comparing best-fit models (Model 2 and 1, respectively, in Table 1) through an Akaike or Bayesian Information Criterion (AIC, or BIC), which penalise models with a larger number of free parameters, clearly favors the model with the more complex vertical structure. Formally, the large BIC differential (or analogously the ratio of the model likelihoods) indicates that model with a double Gaussian vertical structure is $>10^{8}$ times more likely to be better at describing the data than the model with a single Gaussian, providing strong evidence in support of our conclusion.

The narrow Gaussian component has a vertical FWHM of $\sim$2 au at 150 au, which is similar to the size of our resolution element and only resolved through the forward modelling presented here. Note also that the inclination of the double Gaussian model is significantly higher (closer to edge-on) than found for the single Gaussian model. This is likely attributable to the single Gaussian model being driven to a less edge-on configuration in an attempt to capture the broad sky-projected vertical distribution of the disk emission. The position angle East of North of the mm dust disk ($29.7_{-0.1}^{+0.1}$) is marginally consistent within the uncertainties with the outer disk observed in scattered light ($29.3_{-0.1}^{+0.2}$, Lagrange et al., 2012).

4.4 Adding model complexity:

radially varying aspect ratio

The appearance of the residuals in Fig. 6, after subtracting a model where the aspect ratio is the same at all belt radii, indicates that we have already achieved a good fit to the data. Nonetheless, we test whether we can detect a radial dependence of the vertical aspect ratio through a model where the scale height is a function of radius, and this dependence is parametrized through a flaring parameter $\beta$. The dust density distribution reads the same as in Eq. 2, but with $h$ and $h_{2}$ being a function of radius following $h(r)=h_{\rm 50au}\left(\frac{r}{\rm 50au}\right)^{\beta-1}$ and similarly $h_{2}(r)=h_{2,\rm 50au}\left(\frac{r}{\rm 50au}\right)^{\beta-1}$. Then, $\beta=1$ indicates a constant aspect ratio $h=H(r)/r$ as considered in the previous models, $\beta>1$ a flared disk, and $\beta=0$ a disk with a radially constant scale height.

The marginalized posterior probability distribution of $\beta$ results in $\beta=0.7^{+0.2}_{-0.2}$, indicating that models with a scale height increasing with radius ($\beta>0$) are significantly preferred, with a marginal improvement achieved for models with an aspect ratio decreasing with radius ($\beta<1$). However, we find no evidence that the extra free parameter $\beta$ leads to an overall improvement in the fit to the data, which is confirmed by visual inspection of the residuals. We therefore conclude that our data supports models where the scale height increases with radius, with a flaring parameter $\beta$ consistent with one.

4.5 Functional forms other than Gaussians

Our results do not rule out that a different parametric prescription for the vertical density distribution could describe the data as well as the double Gaussian parametrization adopted here. For example, we here attempt to use a model where the dropoff of the vertical density distribution differs from a Gaussian (as implemented by Ahmic et al., 2009). The functional form reads

[TABLE]

with $h(r)=h_{\rm 50au}\left(\frac{r}{\rm 50au}\right)^{\beta-1}$ (as in §4.4), and $C=\left[\int^{+\infty}_{-\infty}{e^{-\left(\frac{\lvert z\rvert}{\sqrt{2}hr}\right)^{p}}}\right]^{-1}$ ensuring normalization of the vertical distribution. Distributions with $p<2$ allow the vertical distribution to drop faster than a Gaussian the further we move from the midplane, but to have broader wings; Ahmic et al. (2009) found this scenario to be preferred in scattered light observations of the $\beta$ Pic belt, with $p$ values around or below 1.

We tried two model runs with the above functional form and $p$ left as a free parameter. In the first (Model 4), we fix the aspect ratio ($\beta=1$) as done for the Model 1 and 2 runs. In the second run (Model 5), we also leave the flaring parameter $\beta$ free. Our fitting results for both model runs, shown in Table 1, indicate that $p\sim 0.9$, in broad agreement with the scattered light findings, and consistent with the broad-winged distribution derived from the image analysis.

We also find that best-fit Model 4 and 5 (with $p=0.9$) are significantly better than a single Gaussian (Model 1). Compared with the double Gaussian, constant aspect ratio model (Model 2, see Table 1, bottom), we find that Model 4 (with $\beta$ fixed to 1) is slightly worse, whereas Model 5 (with $\beta\sim 0.4$) is just as good at reproducing the data, with extremely similar residuals. In summary, there is no evidence in the mm data that a model with a single population of mm grains that has a vertical dropoff rate that is shallower than Gaussian (Eq. 3 with $p\sim 0.9$, Models 4 and 5) is preferred to a double Gaussian model (Model 2). In fact, the convergence of Models 4 and 5 to a broad-winged distribution confirms our main finding that there exists a population of high inclination (dynamically hot) particles in the $\beta$ Pic disk. For later interpretation, we prefer to use the double Gaussian simply to facilitate direct comparison with the Kuiper belt (see §5.2.3).

4.6 On the SW/NE asymmetry in the mm dust: comparison with the 885 $\mu$m data

Independent of the models explored above, the vertically integrated radial profile of the residuals of the 1.33 mm data shows no significant SW/NE asymmetry (bottom left panel of Fig. 5 and 6, where the shaded region is the $\pm 1\sigma$ confidence interval). Averaging emission between 30 and 80 au on each side of the star, as in §3, we find a negligible residual flux difference between the NE and the SW of $82\pm 85$ $\mu$Jy. Varying the ranges over which we radially average does not affect this result. We conclude that we do not detect a significant SW/NE asymmetry in the new 1.33 mm dataset.

Given the discrepancy with the 15% asymmetry reported in the 885 $\mu$m dataset (Dent et al., 2014), we independently re-reduced the dataset as described in §2.2 of Matrà et al. (2017). Just like the 1.33 mm, compact configuration dataset, this consisted of a mosaic pointed $\pm 5\arcsec$ from the stellar location, and along the belt’s PA. We imaged the continuum using the tclean task (CASA v5.1.0) in multi-frequency synthesis, mosaic mode, using natural weights and multiscale deconvolution to best recover faint, extended emission. We obtained an image with synthesized beam size of $0.70\arcsec\times 0.55\arcsec$ at a PA of $64\fdg 4$, and measure an RMS noise level of 70$\mu$Jy beam*-1*. The spatially integrated line flux measured on the primary-beam-corrected map as in the 1.33 mm data is $66\pm 6$ mJy.

As done for the 1.33 mm data, we extract the subset of the dataset containing visibilities at $u-v$ distances $>300k\lambda$ to separate the stellar emission from the disk contribution. We find that the star is not significantly detected. This allows us to set a 3$\sigma$ upper limit of $<786$ $\mu$Jy on the stellar flux at 885 $\mu$m, which is consistent with the expected 196 $\mu$Jy obtained through Rayleigh-Jeans extrapolation from near-IR wavelengths. We therefore include the star with a fixed flux of 196 $\mu$Jy in the modelling below.

We simultaneously fit the visibilities obtained for both pointings as described in §4.1, using the single and double vertical Gaussian models described in §4.2 and §4.3. For each of the two models, we fix the model parameters to the best fit values found for the higher resolution and sensitivity, 1.33 mm dataset (Table 1, Model 1 and 2). We leave as free parameters the scale factor $\rho_{0}$ governing the total flux of the belt, which is different at the shorter wavelength, as well as an RA and Dec offset and a weight-rescaling factor for each of the two pointings that should be different from the 1.33 mm datasets.

Fig. 7 and 8 show residual maps for the single and double vertical Gaussian models, respectively. We find the residuals, particularly the radial and vertical profiles, to be remarkably similar to those obtained for the 1.33 mm observations. In fact, despite the lower resolution and sensitivity, the 885 $\mu$m dataset is also better fitted by a model with a double over a single Gaussian vertical distribution, and also shows no significant SW-NE asymmetry. Averaging emission between 30 and 80 au on each side of the star, as in Dent et al. (2014), we find a negligible residual flux difference between the NE and the SW sides of $60\pm 130$ $\mu$Jy.

An interesting hypothesis is that the two Gaussian populations have different spectral indices, which could be caused by a difference in the size distribution of grains between the two populations. While this can in principle be tested with resolved observations at different wavelengths, such as presented here, the coarse resolution of the 885 $\mu$m data (which, at 100 au from the star, is of order the scale height of the broad Gaussian population) and the relatively small difference in wavelength between the two datasets do not allow us to test this hypothesis. Future ALMA observations at a S/N and resolution comparable to the 1.33 mm dataset presented here, and at a wavelength as far as possible from 1.33 mm, are needed to explore this idea.

5 Discussion

In §3 and §4 we modelled ALMA observations of the $\beta$ Pictoris disk at 1.33 mm and 885 $\mu$m, which show the belt’s vertical dust density distribution is much better fitted by the sum of two Gaussian distribution compared to a single Gaussian distribution. In this Section, we will link this vertical profile to the distribution of particle inclinations, and connect it to their dynamical excitation to infer the potential action of large stirring bodies within or exterior to the belt.

5.1 Connecting a belt’s vertical density distribution to the planetesimals’ inclination distribution

The observed inclination distributions of planetesimals in the Kuiper belt was first investigated by Brown (2001, B01 hereafter). The authors used a Gaussian multiplied by a sine function to model the observed inclinations of Kuiper Belt Objects (KBOs), and showed that this functional form naturally arises from N-body simulations from gravitational perturbations, or stirring, within a planetesimal disk with initially zero inclinations. The function reads $f(i)di\propto\sin{i}\ e^{-\frac{i^{2}}{2\sigma^{2}_{i}}}di$ and for small angles corresponds to a Rayleigh distribution, $f(i)\propto\frac{i}{\sigma_{i}^{2}}e^{-\frac{i^{2}}{2\sigma^{2}_{i}}}$. This result is consistent with other dynamical studies also finding that gravitational perturbations of planetesimals in a thin disk lead to Rayleigh distributions of inclinations and eccentricities (e.g. Lissauer & Stewart, 1993; Ida & Makino, 1992). For direct comparison with the Kuiper belt, we here follow the notation of B01 but note that dynamical studies typically express the inclination distribution as a function of the mean of the squares of the inclinations $\langle i^{2}\rangle$, with $\langle i^{2}\rangle=2\sigma_{i}^{2}$ in the expressions quoted above. Fig. 9 (left) shows an example inclination distribution (following B01) for a population of particles with a mean inclination of $\sim 7.9^{\circ}$ (orange points) and shows how this is well matched to a Rayleigh distribution (blue line).

Under the assumption of circular orbits, this inclination distribution of planetesimal orbits $f(i)$ in a belt at a given radius can be linked to the distribution of latitudes $L(\theta)$ (about a planetary system’s ecliptic plane) at which the planetesimals are found (B01), through

[TABLE]

The integral can be simplified through use of the small angle approximation ($\sin{x}\sim x$ for both $i$ and $\theta$) given the small inclination angles considered here. For example, Fig. 9 (centre and right) shows the latitudinal and height distribution derived from an inclination distribution (left) with a mean $\bar{i}$ of 7.9*∘* and an RMS of 8.9*∘*, corresponding to the hot population of the $\beta$ Pic belt. At the height corresponding to 3$\sigma$ of the distribution (16.5 au, implying a latitude of $18.9^{\circ}$), the small angle approximation is accurate to within 2%. Taking the small angle approximation, the integral in Eq. 4 becomes

[TABLE]

and can be solved analytically leading to the expression

[TABLE]

Given the small inclinations $i$ expected will inevitably lead to small latitudes $\theta$, we can again use the small angle approximation to obtain $\cos{\theta}\sim 1$ and $\mathrm{erf}\left(\sqrt{\frac{\frac{\pi^{2}}{4}-\theta^{2}}{2\sigma_{i}^{2}}}\right)\sim 1$ (the latter easily shown given e.g. $\theta\ll 5\sigma_{i}$ and small $\theta$). Then, we come to the conclusion that the latitudinal distribution of particles with a Rayleigh distribution of inclinations (Fig. 9, left) is simply a Gaussian with standard deviation equal to the $\sigma_{i}$ parameter of the inclination distribution, i.e. $L(\theta)\propto e^{-\frac{\theta^{2}}{2\sigma_{i}^{2}}}$. We verify this in Fig. 9 (centre) by comparing an evaluation of the numerical integral in Eq. 4 with the Gaussian derived analytically above.

Finally, for small angles the distribution of heights (rather than latitudes) of the planetesimals will also be very well approximated by the Gaussian

[TABLE]

where $\sigma_{z}=H=hr\sim\sigma_{i}r$. Thus, using a Gaussian vertical distribution in our model fits of §4 is justified by the expectation of a Rayleigh distribution of inclinations for a belt whose dynamics is dominated by gravitational stirring from bodies within the belt alone (e.g. Lissauer & Stewart, 1993; Stewart & Ida, 2000). The aspect ratio of an observed vertical Gaussian distribution can then be used to derive the mean $\bar{i}$ and root mean square $\mathrm{RMS}_{i}=\sqrt{\langle i^{2}\rangle}$ of the inclination distribution of belt particles, following

[TABLE]

In §3 and 4 we showed how the vertical distribution of mm dust in the $\beta$ Pictoris belt is best described by a sum of two Gaussian distributions. The derivation above then indicates that the corresponding distribution of particle inclinations is bimodal, corresponding to the sum of two Rayleigh distributions (Fig. 10). From now on, we will refer to these as the hot and cold populations of particles, with high and low inclinations respectively. The best-fit RMS inclination derived from Model 2 is $0.156^{+0.011}_{-0.009}$ rad ($8.9^{+0.7}_{-0.5}$ degrees) for the hot population and 0.02${}^{+0.01}_{-0.01}$ rad ($1.1^{+0.5}_{-0.5}$ degrees) for the cold population.

5.2 On the origin of the inclination distribution

5.2.1 External stirring: $\beta$ Pic b

The only known body which could contribute to structure in the belt is the directly imaged, super-Jupiter mass planet $\beta$ Pic b. The most recent orbital determination for the planet indicates a best-fit semimajor axis of 8.8 au, an inclination to the line of sight $I_{\rm p}$ of 89.3*∘* and a position angle $\Omega_{\rm p}$ of 32.1*∘* (Lagrange et al., 2018a). When comparing to the belt as resolved by ALMA, we confirm that the orbital plane of the planet and the belt are misaligned with respect to one another. This is thanks to our high resolution measurement of the belt’s position angle $\Omega_{\rm b}$ and inclination to the line of sight $I_{\rm b}$ from thermal emission, which for the first time are free of assumptions on the unknown phase function of the grains, which affected previous determination from scattered light observations.

The misalignment of the orbital plane of the planet and that of the disk can be characterized by the angle $\mathrm{d}I$, the inclination of the planet orbit with respect to the disk midplane. A simple transformation between reference frames shows that

[TABLE]

We find that the belt’s inclination $I_{\rm b}$ to the line of sight (see Table 1) for all the best-fit models to be consistent with that of $\beta$ Pic b’s orbit ($I_{\rm p}\sim$89.3 deg, but see Fig. 5 in Lagrange et al., 2018a, for an estimate of the uncertainty). If confirmed by future observations, that would indicate a planet-belt orbital plane misalignment $\mathrm{d}I\sim 2.4$ deg, approximately equal to the difference in position angle between the belt and the planet’s orbital plane.

Assuming this misalignment has been retained for the $\sim$23 Myr age of the system, the misaligned orbital plane of $\beta$ Pic b should significantly affect the belt’s vertical structure through secular perturbations affecting the belt’s planetesimals. In particular, $\beta$ Pic b should impose a forced inclination $i_{\rm p}=\mathrm{d}I\sim 2.4$ deg on the belt’s particles, and cause their inclinations to oscillate between 0 and 2$i_{\rm p}$ (e.g. Dawson et al., 2011; Nesvold & Kuchner, 2015). Particles further out in the belt are affected on longer timescales; this has long been considered the origin of the belt’s warp observed in the inner regions out to $\sim$100 au from scattered light observations (Mouillet et al., 1997; Augereau et al., 2001), that trace the smallest grains in the collisional cascade.

In this scenario, we would expect 1) the warp to be also present for mm grains, and 2) the vertical distribution of dust density, and hence mm emission, to be enhanced at 0 and 2$i_{\rm p}$ where particles spend most of their time during their oscillations. An example of the expected disk morphology imparted by $\beta$ Pic b is shown in Fig. 1 of Nesvold & Kuchner (2015). This shows clear vertical flux enhancements displaced $\sim\pm 9$ au from the midplane at 100 au from the star, which are clearly ruled out by the disk spine and vertical profiles derived from the new ALMA data (Fig. 3, central and top panel). Indeed, the mm continuum emission is centrally peaked and displaced by no more than 5 au vertically (at the $3\sigma$ level) throughout the disk midplane. This would limit any sky-projected warp, if present, to less than 2.8 deg (when measured at 100 au), which roughly corresponds to the forced inclination imposed by the planet.

In summary, we find that the disk spine and centrally peaked vertical emission profile are inconsistent with models of the secular perturbations imposed by $\beta$ Pic b on the belt. Taken alone, this could indicate that either 1) the planet is less inclined compared to the belt than currently believed, 2) another massive planet is present in the outer regions of the belt, dominating its dynamics, or 3) $\beta$ Pic b has only recently been put on a misaligned orbit, and has not had time to interact with the outer belt. Such an alternative picture would have to be reconciled with the scattered light and CO inner disk tilt, and with the hot and cold inclination populations reported here.

5.2.2 Gravitational stirring within the belt and constraints on the size of the belt’s largest bodies

The strongest constraint on the double-Gaussian vertical distribution of the belt comes from its radial outer edge, which is least affected by line-of-sight integration of emission originating at different orbital radii. Even if the belt is misaligned with $\beta$ Pic b’s orbit, this outer region has not had enough time to be affected by secular perturbations from the planet within the age of the system. Therefore, any structure in this region, and particularly the hot population, is either inherited from the gas-rich protoplanetary phase of evolution (e.g. Walmswell et al., 2013), or must be the result of dynamical excitation from a large body/bodies other than $\beta$ Pic b.

We here consider gravitational stirring by large bodies within the belt itself. Encounters with such bodies will act as a source of dynamical heating, setting the velocity dispersion of particles, and producing a single, Rayleigh distribution of inclinations and eccentricities (see Kokubo & Ida, 2012, and references therein). Given that the observed inclination distribution of the $\beta$ Pic belt clearly deviates from a single Rayleigh distribution, we can rule out gravitational stirring from within the belt alone as the source of the observed vertical structure. Nonetheless, we here consider gravitational stirring as a potential source of either the hot or cold populations, while keeping in mind that other mechanisms are required to produce the other population.

Gravitational stirring leads to a relative velocity of particles of the form $v_{\rm rel}=\sqrt{1.25\langle e^{2}\rangle+\langle i^{2}\rangle}v_{\rm Kep}$, where $v_{\rm Kep}$ is the local Keplerian velocity (e.g. Lissauer & Stewart, 1993). Since $\sqrt{\langle e^{2}\rangle}=2\sqrt{\langle i^{2}\rangle}$ (Ida & Makino, 1992), the RMS of the inclination distribution is a direct measure of the relative velocity of particles, following

[TABLE]

where $h$ is unitless, $r$ is in au, $M_{\star}$ is in Solar masses, and the resulting $v_{\rm rel}$ is in km/s. Then, for an assumed stellar mass of 1.75 M*⊙* (Crifo et al., 1997), the cold and hot populations in the $\beta$ Pic belt (Model 2) indicate approximate relative velocities of 0.27 and 2.1 km/s at 50 au, dropping to 0.16 and 1.23 km/s at 150 au.

Viscous stirring from a large body (or bodies) will act to increase the velocity dispersion of planetesimals, and consequently their inclinations, over time. Following Ida & Makino (1993), assuming the system has been evolving from a perfectly cold disk ($e,i\sim 0$ at $t=0$ Myr) for the age of the system ($t\sim 23$ Myr) we can then relate the measured inclinations to the product of the mass and the surface density of the largest bodies.

We assume we are in the dispersion-dominated regime and that stirring is dominated by large body/bodies rather than the small planetesimals themselves. Combining Eq. 4.2 and 4.9 of Ida & Makino (1993) and solving for the inclination as a function of time, we obtain

[TABLE]

where the constant $C_{\rm I}\sim 2$ (Ida & Makino, 1993), $\Omega$ is the Keplerian frequency in rad s*-1*, $M_{\star}$ is the stellar mass in kg, $t$ is time in s and $\Sigma$ is the surface mass density of the large stirring bodies (in kg m*-2*), each of mass $M$ (in kg). The measured RMS inclination of a planetesimal population then allows us to set a constraint on the mass times the surface density of the largest bodies of

[TABLE]

(where $M_{\star}$ is now in M*⊙, $r$ is in au, and $t$ is in Myr). This leads to a value of $5.9\times 10^{-6}$ M${}_{\oplus}^{2}$au-2* for the cold population, and 0.022 M${}_{\oplus}^{2}$au*-2* for the hot population, at 150 au in the $\beta$ Pic belt. This corresponds to the purple lines in the $\Sigma$ vs $M$ space shown in Fig. 11.

Producing the observed excitation however requires the bodies causing the stirring to not collide (and hence be destroyed) within the age of the system. The rate of mutual collisions for a population of bodies of equal mass $M$ can be estimated as

[TABLE]

where $\sigma=4\pi R^{2}$ is the effective cross section of the collision (with $R=(0.24M/\rho)^{1/3}$ being the radius of the large bodies, and $\rho$ their material density), $n$ is the number density of such bodies, and $v_{\rm esc}=\sqrt{GM_{\star}/R}$ their escape velocity. As in Eq. 8, we assume $v_{\rm rel}=\sqrt{6\langle i_{p}^{2}\rangle}$ where $\langle i_{p}^{2}\rangle$ is the RMS inclination of the large bodies doing the stirring, which is likely primordial and here assumed to be 0.05. Setting the collision timescale ($t_{\rm c}=R_{\rm col}^{-1}$) to be greater than the age of the system, we can obtain an upper limit to the surface density of the large bodies given their mass, shown as the orange downward triangles in Fig. 11. In log-log space, the slope of this surface density upper limit transitions from positive at the low mass end where $v_{\rm esc}\ll v_{\rm rel}$, to negative at the high mass end where $v_{\rm esc}\gg v_{\rm rel}$.

The intersection between the orange and purple lines indicate the minimum mass necessary for large bodies to produce the observed level of stirring, but at the same time avoid potentially destructive collisions in 23 Myr. It can be shown analytically by combining Eq. 9 and 11 at the high mass end (where $v_{\rm esc}\gg v_{\rm rel}$) that requiring the collision timescale to be equal to the viscous stirring timescale leads to the planetesimals’ relative velocities to be not exactly equal to, but of order the escape velocity of the largest bodies (blue dash-dotted line). This means that this minimum mass can be approximately estimated by simply setting the planetesimals’ relative velocities to be equal to the stirring body’s escape velocity.

A further requirement for the surface density of large bodies in the belt is that at least one large body is present at the radii where the stirred population is observed. In other words, if we consider the stirred population near, say, 150 au, there needs to be at least one large body such that its stirring zone encompasses 150 au. This sets a lower limit on the mass surface density (green upward triangles in Fig. 11),

[TABLE]

where $\Delta r_{\rm stir}=8\sqrt{3}r[M/(3M_{\star})]^{\frac{1}{3}}$ is the width of the stirring zone (equal to twice the large body’s feeding zone, following Eq. 3.4 of Ida & Makino (1993)). The intersection of the purple line and green triangles represents the maximum mass for a body to stir a given population and for at least one such body to exist so that its stirred region encompasses a given radius at which the population is observed.

Finally, we can consider the lack of observed radial gaps in the belt as a further constraint on the large body’s mass. This is because a body sufficiently large will open a gap of width roughly equal to its chaotic zone ($\Delta r_{\rm chaos}\sim 3r(M/M_{\star})^{\frac{2}{7}}$, e.g. Wisdom, 1980; Morrison & Malhotra, 2015; Marino et al., 2018). Simply setting the width of the chaotic zone to be less than the resolution of our observations leads to an upper limit to the mass of the bodies within the belt. At 150 au, we find this mass to be $<$0.15 M*⊕. However, we need to consider the timescale necessary to open this gap. Following Eq. 7 and 8 in Morrison & Malhotra (2015), the timescale for a 0.15 M⊕* body to open a gap at 150 au in the $\beta$ Pic belt is $\sim 5.7$ Gyr. Then, the lack of radial gaps can not only be caused by the chaotic zone not being resolvable, but also by the large body’s gap-opening timescale being longer than the age of the system, where the latter sets a more conservative constraint on the maximum planet mass. At 150 au (50 au), we find that a stirring planet should be less massive than 75 (25) M*⊕* to not produce observable gaps in 23 Myr (brown left-pointing triangles in Fig. 11).

In reality, this estimate is uncertain because the clearing timescale of Morrison & Malhotra (2015) is defined as the ‘time required to reach 50% of the final survivor fraction within the cleared zone’ and is therefore not to the timescale to produce an observable gap. The latter may be shorter, for example, if the observing sensitivity is sufficient to detect gap-belt surface brightness contrasts corresponding to a higher survival fraction of particles in the belt. At the same time, the edge-on viewing geometry would make it harder to detect such gap-belt contrast compared to a face-on belt, potentially making the timescale to produce a detectable gap longer. Detailed N-body simulations combined with radiative transfer and observing simulations are needed to accurately pinpoint this timescale (see e.g. Marino et al., 2018), but are beyond the scope of this work.

Combining the dynamical constraints, we find that the mass of the largest bodies within the belt necessary to stir a hot population such as observed for $\beta$ Pic over 23 Myr should be between 0.007 and 54 M*⊕* at 150 au. On the other hand, if we consider viscous stirring as the origin of the cold population, we find that the largest bodies must have masses below 0.4 M*⊕* at 150 au (or 0.15 M*⊕* at 50 au) to have maintained such a low dynamical excitation. These masses are all far below and therefore consistent with current detection limits from direct imaging (Lagrange et al., 2018b).

To conclude, we underline that viscous stirring by large bodies within the belt alone cannot produce the observed complexity of the vertical structure of the $\beta$ Pictoris belt. However, it could produce either the hot or cold population, while requiring another mechanism to produce the other. If that were to be the case, the observed dynamical excitation of the hot (cold) populations could be produced through viscous stirring by bodies of mass $0.007-54$ ($<0.4$) M*⊕* at 150 au.

5.2.3 Outward migration of a $\beta$ Pic c

A bimodal distribution of inclinations such as observed in the $\beta$ Pictoris belt was discovered for our own Kuiper belt (KB) by B01 (after debiasing was taken into account, see also e.g. Kavelaars et al., 2009; Petit et al., 2011). For both belts, we clarify that a single, broad non-Rayleigh distribution of inclinations cannot be excluded (e.g. Volk & Malhotra, 2011); we just follow the original, bimodal functional form of B01 for direct comparison.

Fig. 12 shows how the inclination distribution of the $\beta$ Pic belt (namely, best-fit Model 2) compares with that of the KB from B01. Overall, we find lower inclinations for both the hot and cold inclination populations of the $\beta$ Pic belt when compared to the KB, indicating an overall lower level of dynamical excitation. The fraction of all particles residing in the cold population are however consistent between the $\beta$ Pic belt and the KB ($20^{+5}_{-4}$% and $\sim 19-26$%), with both belts’ mass being dominated by the hot population.

In the Kuiper belt, the high inclinations seen for the hot classicals as well as the resonant and scattered populations set important constraints on its formation scenario (for a review, see e.g. Nesvorný, 2015). These populations likely formed in a massive planetesimal disk within 30 au and then got implanted into their present day orbit by outward-migrating Neptune. The high inclinations can be attained through sweeping resonances (e.g. Levison & Morbidelli, 2003) or scattering by Neptune during its migration (Gomes, 2003). Nesvorný (2015) showed that the latter mechanism can work as long as Neptune’s migration was slow ($>$10 Myr).

By analogy with the Kuiper belt, we consider whether $\beta$ Pic’s structure could be produced by outward migration of a hypothetical $\beta$ Pic c planet located at the inner edge of the belt. The presence of such a planet has been suggested to produce the SW/NE asymmetry observed for the CO gas (Dent et al., 2014; Matrà et al., 2017) as well as for the small dust grains at wavelengths below 24 $\mu$m (e.g. Telesco et al., 2005), following models predicting resonant planetesimals to concentrate at specific longitudes relative to the planet (Wyatt, 2003). To produce the asymmetry observed for the CO, two conditions need to be met. First, the 2:1(u) resonance trapping probability (where the (u) indicates planetesimals in 2:1 resonance whose libration center increases from 180*∘* during the migration) needs to be sufficiently high, and second, the planetesimal eccentricities need to be excited to a high enough level (we here require e$>$0.3, see Fig. 6, top row of Wyatt, 2003).

These eccentricities can be excited by resonant forces, causing them to increase as the radial extent of the migration increases (Eq. 22 in Wyatt, 2003). We find that the ratio between the semimajor axis of planetesimals after and before migration should be $\gtrsim 1.8$ for eccentricities to reach values above 0.3. For example, outward migration of a planet from 31.4 to at least 57.4 au (close to the inner edge of the belt as observed today) could have swept planetesimals in the 2:1(u) resonance from 50 to at least 91 au and produced sufficiently high eccentricities that an asymmetry may be observed. Given this migration must have taken place for at most the 23 Myr age of the system, we estimate a lower limit on the migration rate of 1.13 au/Myr.

Ensuring the 2:1(u) trapping probability is sufficiently high (we here require it to be above 50%, see Eq. 8 and Table 1 of Wyatt, 2003) allows us to link this migration rate to the mass of the migrating planet. We find that a planet with mass of at least $28$ M*⊕* must have migrated at least $\sim$25 au to produce an asymmetry as observed for the CO. Therefore, it is possible that outward migration of a super-Neptune sized planet may have produced the SW/NE asymmetry observed for the CO and for the small dust at short wavelengths. Can this however be reconciled with the lack of asymmetry for the mm-sized dust?

Wyatt (2006) showed that the appearance of resonant structure depends on the size of the observed grains, with grains below a critical size $D_{\rm crit}$ falling out of resonance. This size depends solely on the migrating planet’s mass (and some stellar parameters, see Eq. 14 in Wyatt, 2003). For a planet of at least $28$ M*⊕*, this size corresponds to 242$\mu$m, although grains up to $\sim$10 times this size can show significantly smeared resonant structure. The observed emission in the ALMA 1.33mm dataset is dominated by grains of size similar to the observing wavelength or a few times smaller (see e.g. Fig. 6, bottom left in Wyatt, 2003). It therefore seems plausible that the resonant structure could be smeared out to the extent that it is undetected in the 1.33mm ALMA data as well as in the 24$\mu$m (Telesco et al., 2005) and 70-500 $\mu$m (Vandenbussche et al., 2010) datasets. Observations at wavelengths even longer than presented here could clarify whether this is the origin of the lack of observed structure at 1.33mm.

The model predicts the structure of the short lived, smallest grains and the CO to be asymmetric. This means that the asymmetric resonance structure is only lost, or smeared out, at intermediate grain sizes (corresponding to mid to far-IR/mm wavelengths). Asymmetric resonant structure should therefore be seen for the smallest grains, for the CO, and for the large resonant grains of sizes $\gg 242\mu$m. Note, however, that while the distribution of CO and blow-out grains traces the collisional mass loss rate from the belt (which is proportional to the relative velocity of particles multiplied by the square of their cross section per unit volume), the distribution of the large resonant grains simply traces the particles’ cross sections per unit volume. That has the important implication that the SW-NE asymmetry will be much less pronounced and potentially hardly detectable for large resonant grains when compared to the blow-out grains and the CO. This means that the ALMA 1.33 mm and any longer wavelength data could probe resonant grains, but the resonant structure may be too faint for detection.

As well as the SW/NE asymmetry and the high inclinations observed at the estimated location of the resonant clumps, the outward-migrating planet needs to be able to produce high inclinations all the way out to the outer edge of the disk (i.e. beyond the resonant clumps), as well as the tilt observed for the CO and blow-out grains. A hot population far out in the disk may be produced via outward scattering, as seen for the scattered population of the Kuiper Belt (Brown, 2001). At the same time, the inner tilt seen for CO and small grains may be a projection effect caused by two asymmetric, azimuthal resonance clumps in front and behind the plane of the sky in a non perfectly edge-on configuration (see e.g. Matrà et al., 2017, Fig. 8, bottom).

Overall, this $\beta$ Pic c migration scenario could qualitatively reproduce the observed features of the outer belt, but needs to be tested through future planet migration simulations. These should be tailored to reproduce $\beta$ Pic’s gas and dust structure as seen by ALMA, and should take into account the influence that $\beta$ Pic b would have on such a planet, which we did not consider here.

6 Conclusions

In this work, we presented new ALMA observations of the $\beta$ Pictoris belt at mm wavelengths. This provided us with the highest resolution picture of large, bound grains in the belt to date, allowing us to study their vertical distribution for the first time. We report the following findings, confirmed by both image analysis and modelling of the interferometric visibilities from the ALMA data:

• •

The vertical distribution of mm grains cannot be fitted by a single Gaussian. In turn, this implies that a single Rayleigh distribution is not a good description of the particles’ inclination distribution. We find a model with a bimodal distribution of inclinations, analogous to the Kuiper belt, to be a much better fit to the data. RMS inclinations are $0.156^{+0.011}_{-0.009}$ rad (8.9${}^{+0.7}_{-0.5}$ degrees) for the hot component, and $0.02^{+0.01}_{-0.01}$ rad (1.1${}^{+0.5}_{-0.5}$ degrees) for the cold component, with the hot component containing $80^{+4}_{-5}$% of the observed mass.

• •

The position angle and inclination of the belt to the line of sight, here determined free of bias from the scattered light phase function, indicate significant misalignment between the belt’s and $\beta$ Pic b’s orbital planes. Focusing on the belt’s inclination to the line of sight, we find that it is consistent with that of $\beta$ Pic b’s orbit, which if confirmed would indicate that misalignment is maximum in or close to the plane of the sky.

• •

After fitting an axisymmetric model to the new 1.3 mm and archival 885 $\mu$m visibilities, we find no significant evidence for the previously reported SW/NE brightness asymmetry in the continuum emission arising from mm grains. This is in agreement with previous 24-500$\mu$m dust observations, but in stark contrast with the asymmetries seen in both gas and optical/near-infrared observations.

• •

Similarly, we find no significant evidence for a tilt of the inner disk (i.e. a warp) as pronounced as observed in scattered light and CO observations. While the presence of a smaller tilt is not ruled out, this would now require a mechanism to produce a larger tilt in scattered light and CO gas with respect to the mm grains.

We then consider possible origins for the observed inclination distribution observed in the mm grains. We begin by considering the misalignment of $\beta$ Pic b with respect to the belt, which should force the particles’ inclination to oscillate between 0 and twice the planet’s inclination with respect to the belt plane. This is ruled out by the lack of significant warping and by the vertical distribution observed by ALMA. Furthermore, the strongest constraint on the presence of the hot population comes from the belt’s outer edge, which cannot be affected by $\beta$ Pic b’s secular perturbations within the age of the system.

We therefore explore alternative scenarios, such as viscous stirring from large bodies within the belt. The latter alone cannot explain the observed inclination distribution, but could be partly responsible for the excitation of either the cold or hot populations. We lay down a framework to constrain the mass and surface density of large bodies by connecting the observed inclination distribution to the expectation from viscous stirring. We find that producing the hot population would require large bodies of mass between 0.007 and $\sim$54 M*⊕* at 150 au within the belt. Most of these bodies cannot produce observable gaps as their gap-opening timescale is longer than the age of the system, although whether the most massive ones $\gtrsim 10$ M*⊕* could is sensitive to the adopted ’gap opening’ definition. On the other hand, keeping the low dynamical excitation of the cold population would require bodies to have masses below 0.4 M*⊕* at 150 au.

Another potential explanation could come from the presence of an outward-migrating $\beta$ Pic c at the inner edge of the belt. Such a planet could produce both the SW/NE asymmetry seen in CO and small grains through resonant sweeping, and explain the hot population of inclinations observed in a way analogous to the inferred outward migration of Neptune in the Solar System. A $\gtrsim 28$ M*⊕* planet having migrated from $\sim$31 to 57 au in the system could produce the SW/NE asymmetry for the CO and small grains. This can likely be reconciled with the lack of a SW/NE asymmetry for grains observed between 24$\mu$m and 1.33 mm, but requires tailored dynamical simulations, for example accounting for the influence of $\beta$ Pic b on such a planet, to draw a more robust conclusion.

To conclude, new evidence from ALMA data combined with dynamical arguments indicate that $\beta$ Pic b is unlikely to be the only large body responsible for the belt’s observed vertical structure, unless planetesimals were born on high inclination orbits in the earlier protoplanetary disk phase of evolution. Instead, we have shown that the presence of a hot population of inclinations, analogous to the Kuiper belt’s, could be produced through dynamical excitation by additional large bodies interior and/or exterior to the belt. Detailed dynamical simulations combined with observations at longer wavelengths or higher resolution are necessary to confirm or rule out the scenarios proposed here.

LM acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. GMK is supported by the Royal Society as a Royal Society University Research Fellow. This paper makes use of ALMA data ADS/JAO.ALMA#2012.1.00142.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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. ALMA.