ALMA observations of the multiplanet system 61 Vir: What lies outside super-Earth systems?

TL;DR

This study uses ALMA and JCMT data to resolve the debris disc around 61 Vir, revealing its extent, surface density profile, and implications for unseen planets influencing its structure.

Contribution

First resolved millimetre imaging of 61 Vir's debris disc, combined with a collisional model, to infer disc properties and potential unseen planetary perturbers.

Findings

Disc extends from 30 to at least 150 AU.

Surface density profile has a shallow power law slope.

Unseen planet likely more massive than 10 Earth masses between 10-20 AU.

Abstract

A decade of surveys has hinted at a possible higher occurrence rate of debris discs in systems hosting low mass planets. This could be due to common favourable forming conditions for rocky planets close in and planetesimals at large radii. In this paper we present the first resolved millimetre study of the debris disc in the 4.6 Gyr old multiplanet system 61 Vir, combining ALMA and JCMT data at 0.86 mm. We fit the data using a parametric disc model, finding that the disc of planetesimals extends from 30 AU to at least 150 AU, with a surface density distribution of millimetre sized grains with a power law slope of 0.1$^{+1.1}_{-0.8}$. We also present a numerical collisional model that can predict the evolution of the surface density of millimetre grains for a given primordial disc, finding that it does not necessarily have the same radial profile as the total mass surface density (as previous studies suggested for the optical depth), with the former being flatter. Finally, we find that if the planetesimal disc was stirred at 150 AU by an additional unseen planet, that planet should be more massive than 10 M$_{\oplus}$ and lie between 10-20 AU. Lower planet masses and semi-major axes down to 4 AU are possible for eccentricities $\gg$ 0.1.

Tables

Table 1. Table 1: Best fit values of the ALMA and SCUBA2 data combined. Median ± plus-or-minus \pm uncertainty based on the 16th and 84th percentile of the marginalised distributions.

Parameter	Best fit value
$R_{\max }$ [AU]	${153}_{-26}^{+39}$
$\alpha$	${1.2}_{-1.4}^{+2.3}$
$F_{\mathrm{disc}}$ [mJy]	${3.7}_{-1.1}^{+1.2}$
PA [∘]	$59\pm 5$
$i$ [∘]	$82\pm 4$

Equations

$${\chi }_{\mathrm{ALMA}}^2=\sum _i\frac{||V_{\mathrm{data},\mathrm{i}}-V_{\mathrm{model},\mathrm{i}}|{|}^2}{\delta V_{\mathrm{data},\mathrm{i}}^2},$$

$$Q_{\mathrm{D}}^{\star }=\left[Q_{\mathrm{D},\mathrm{s}}{\left(\frac{\mathrm{D}}{1\mathrm{m}}\right)}^{b_{\mathrm{s}}}+Q_{\mathrm{D},\mathrm{g}}{\left(\frac{\mathrm{D}}{1\mathrm{m}}\right)}^{b_{\mathrm{g}}}\right]{\left(\frac{v_{\mathrm{rel}}}{v_0}\right)}^{1/2},$$

$$M(r,t,D)=M_0^{\prime }(D)\frac{t_{\mathrm{c}}(0)}{t_{\mathrm{age}}},$$

$$q_i=\frac{21+b_i}{6+b_i},$$

$${\Sigma }_{\mathrm{mm}}(r)\propto {\left[r^2{\Sigma }_0(r)\right]}^{\frac{2+k_2-k_2{q}_2}{2+q_2-q_3+k_2-k_2{q}_2}}{r}^{-2+\frac{(19+2q_2)(q_2-q_3)}{6(2+q_2-q_3+k_2-k_2{q}_2)}},$$

$${\Sigma }_{\mathrm{mm}}(r)\propto r^{0.6\alpha +0.9}.$$

$$\begin{array}{r}{\Sigma }_{\mathrm{mm}}(r>r_c)\approx 2{\left(\frac{r}{1\mathrm{AU}}\right)}^{0.6\alpha +0.9}{\left(\frac{t}{1\mathrm{Gyr}}\right)}^{-0.4}{\left(\frac{D_{\max }}{100\mathrm{km}}\right)}^{-0.1}\\ {\left(\frac{{\Sigma }_0}{1\mathrm{MMSN}}\right)}^{0.6}{\mathrm{M}}_{\oplus }{\mathrm{AU}}^{-2},\end{array}$$

$$r_c\approx 4{\left(\frac{t}{1\mathrm{Myr}}\right)}^{\frac{1}{-\alpha +1.5}}{\left(\frac{{\Sigma }_0}{1\mathrm{MMSN}}\right)}^{\frac{1}{-\alpha +1.5}}{\left(\frac{D_{\max }}{100\mathrm{km}}\right)}^{-0.5}\mathrm{AU},$$

$$a_{\mathrm{plt}}(1-e)-5R_{\mathrm{H},\mathrm{q}}>a_{\mathrm{c}}(1+e_{\mathrm{c}}),$$

Full text

††pagerange: ALMA observations of the multiplanet system 61 Vir: What lies outside super-Earth systems?–ALMA observations of the multiplanet system 61 Vir: What lies outside super-Earth systems?

1 Introduction

Planetary systems around main sequence stars are not only composed of planets; planetesimal belts can be also present, analogous to the Kuiper belt (at tens of AU) and the Asteroid belt (within a few AU) in the Solar System. These belts can produce dusty debris discs as the result of a so-called collisional cascade (e.g., Dominik & Decin, 2003; Wyatt et al., 2007), where solids in a wide size distribution from $\mu$m-sized grains up to km-sized planetesimals are ground down in collisions, sustaining high levels of dust and infrared excess over Gyr timescales. Debris discs at tens of AU are fairly common around FGK stars, with occurrence rates of at least $\sim 20\%$ (e.g., Su et al., 2006; Hillenbrand et al., 2008; Carpenter et al., 2009; Eiroa et al., 2013; Thureau et al., 2014; Matthews et al., 2014a; Montesinos et al., 2016); therefore, a complete understanding of their properties can give us information about planet formation and migration history of planets in these systems (e.g., Wyatt, 2006).

Moreover, a few systems are known to host both a planet(s) and a debris disc(s). Among the best studied are $\beta$ Pic (e.g., Smith & Terrile, 1984; Lagrange et al., 2009; Dent et al., 2014), HR 8799 (e.g., Marois et al., 2008; Marois et al., 2010; Matthews et al., 2014b; Booth et al., 2016) and Fomalhaut (e.g., Kalas et al., 2008), all with planets directly imaged and lying between the star and the disc. However, these systems are outliers in terms of their planets and disc properties and neither represents the bulk of the known planetary systems, nor debris discs.

Thanks to unbiased debris disc surveys of FGK stars within 45 pc (e.g., DUNES and DEBRIS Eiroa et al., 2013; Matthews et al., 2014a), it has been possible to study the frequency of circumstellar material around stars hosting high- and low-mass planets detected by radial velocity surveys. Studies focused on high-mass planets found no evidence of a different debris disc incidence rate in these planet hosting stars compared to normal field stars (Greaves et al., 2004; Moro-Martín et al., 2007; Bryden et al., 2009). On the other hand, two studies focused on planetary systems with planet masses below $\sim 95$ M⊕found debris disc incidence rates of: 4/6 (Wyatt et al., 2012), significantly higher compared to field stars; and 2/6 (Moro-Martín et al., 2015), consistent with field stars. Combining both samples, at least four out of eight systems with low-mass planets also have a debris disc, which suggests that there might be a difference in the occurrence of bright debris discs in systems with low-mass planets, as predicted by planet formation models (e.g., Raymond et al., 2011).

One of these planetary systems hosting a debris disc is 61 Vir. This system located at 8.6 pc (van Leeuwen, 2007) is composed of : 1) a G5 $4.6\pm 0.9$ Gyr old star (Wright et al., 2011; Vican, 2012); 2) three RV planets of minimum masses 5, 18 and 23 M⊕ and semi-major axes of 0.05, 0.22 and 0.49 AU, respectively (Vogt et al. 2010, the third one was not confirmed in the HARPS data, Wyatt et al. 2012; and 3) a debris disc discovered by Spitzer with a fractional luminosity $L_{\rm disc}/L_{\star}$ of $2\times 10^{-5}$ (Bryden et al., 2006). The disc was later imaged by Herschel showing that the disc density peaks between 30 and 100 AU and it is inclined by $\sim 77\degr$ with respect to the plane of the sky (Wyatt et al., 2012). If disc and orbits of these planets are co-planar, then the planet masses would be only underestimated by 3%. This system is particularly interesting as the fraction of stars with super-Earths, similar to 61 Vir, could be up to 30-50% (e.g., Howard et al., 2010; Mayor et al., 2011; Fressin et al., 2013), which makes 61 Vir a good case to study the formation of such abundant planets by analysing its debris disc.

Due to a low 50 AU resolution, Herschel could not constrain the exact morphology and dust distribution at the inner regions of the disc, but by image and spectral energy distribution (SED) fitting Wyatt et al. (2012) found three best fit models: 1) an extended disc with a sharp inner edge at $\sim$30 AU, extending at least out to 100 AU, and a surface density or optical depth radial profile with an exponent of -1.1; 2) similar to the first model, but adding an inner component where the surface density increases with radius as $r^{7/3}$ (inspired by collisional evolution models) from 1 AU to the disc inner edge now placed at 43 AU; 3) a two belt model consisting of two 10 AU wide dusty belts centered at 40 and 90 AU. These three models could well fit the previous observations, but the low 50 AU resolution hindered determining the exact dust distribution. Moreover, because the disc emission at Herschel wavelengths is dominated by small grains that are subject to radiation forces, the derived distribution does not necessarily trace the location of the parent planetesimal belt, as they can extend to larger radii beyond that belt (e.g., Thébault & Augereau, 2007).

In this paper we present the first observations of 61 Vir with the Atacama Large Millimeter/submillimeter Array (ALMA) at 0.86 mm, obtained with the aim of studying its debris disc to reveal the location of the parent planetesimals, and place constraints on the presence of planets at large separations that can shape the mass distribution in the disc. Because radiation forces are negligible for mm-sized grains, their distribution can be used to trace the location of the biggest km-sized planetesimals (or bigger), which contain the bulk of the disc mass and sustain the collisional cascade. At millimetre wavelengths, the dust thermal emission is dominated by mm-sized grains ($\sim 0.1-10$ mm), therefore, observations with ALMA are well suited to study the dynamics and origin of debris discs. In order to obtain the best disc constraints, in our analysis we combine new ALMA band 7 observations and new data at 0.85 mm from the Sub-millimetre Common-User Bolometer Array 2 (SCUBA2) installed in the James Clerk Maxwell Telescope (JCMT), thus, incorporating information from small and large angular scale structure.

In addition, we implement a simple numerical collisional evolution model that simulates the evolution of a broad disc, taking into account the disruption threshold of planetesimals as a function of size, how relative velocities vary with radii, and the different features in the size distribution of solids, e.g. the ripples close to the blow-out size. We use this to constrain the initial solid mass or surface density in the disc and the maximum planetesimal size.

This paper is organised as follows. In Sec. 2 we present new SCUBA2/JCMT data. In Sec. 3 we describe the ALMA observations, studying the dust continuum and how it compares with previous Herschel observations. We also search for CO v=0, J=3-2) line emission. In Sec. 4 we fit the SCUBA and ALMA data using a parametric disc model to study the distribution of millimetre dust in the disc. Sec. 5 describes a numerical model to calculate the collisional evolution of a disc at different radii that can be used to compare with observations determining the maximum planetesimal size in a disc and the initial sold mass. In Sec. 6 we discuss the observations and possible scenarios than could explain the low initial solid mass and maximum planetesimal size. Moreover, we constrain the mass, semi-major axis and eccentricity of a hypothetical planet stirring the disc. Finally, In Sec. 7 we summarise and present the main conclusions of this paper.

2 SCUBA2 observations

As part of the SCUBA-2 Observations of Nearby Stars (SONS) survey (Panić et al., 2013), 61 Vir was observed at 0.85 mm with SCUBA2/JCMT (Holland et al., 2013) to constrain the millimetre flux and extent of its debris disc. 61 Vir was observed for 7.5 h and the data was reduced using the Dynamic Iterative Map-Maker within the Starlink SMURF package (Chapin et al., 2013), which was called from the automated pipeline ORAC-DR (Cavanagh et al., 2008). More details on the SCUBA2 data reduction of the SONS survey can be found in Matthews et al. (2015) and Kennedy et al. (2015).

Herschel and VLA observations previously found three background sources close to 61 Vir which could affect our analysis of the SCUBA2 data. To obtain a non-contaminated large scale image and photometry of 61 Vir we subtract these as point sources, using as PSF the SCUBA2 reduced observation of Uranus obtained in the same run. As two of these sources are detected in the ALMA data (see Sec. 3.1), we can derive their fluxes and astrometric positions at 0.86 mm, and accurately subtract these from the SCUBA2 reduced image, correcting for the proper motion of 61 Vir $\mu$=(1.07, -1.06) $\arcsec$ yr*-1* (van Leeuwen, 2007). The third background and more distant source from 61 Vir is not detected with SCUBA2, and lies outside the field of view of the ALMA observations. In Figure 1 we present the SCUBA2 image smoothed with a Gaussian kernel of FWHM $6\aas@@fstack{\prime\prime}5$ after subtracting the two background sources. Integrating all the emission inside a circumference of $15\arcsec$ radius we find a total flux of $5.0\pm 1.2$ mJy (including the stellar emission and calibration uncertainty), slightly lower but consistent within errors with the previous data presented in Panić et al. (2013).

3 ALMA Observations

ALMA band 7 (0.86 mm) observations of 61 Vir were carried out on 2015 April, split into 4 scheduling blocks (one on April 9 and three on April 22) as part of the project 2013.1.00359.S (PI: M.C. Wyatt). The total number of antennas was 44, with baselines ranging from 15 to 349 m, with 5th and 95th percentiles equivalent to 29 and 228 m. This allows us to recover angular scales of $0\aas@@fstack{\prime\prime}6$ up to $6\arcsec$ on the sky.

The correlator was set up with three spectral windows to image the continuum centered at 333.84, 335.78 and 347.74 GHz, each with 128 channels and a total bandwidth of 1.88 GHz; and a fourth one to search for CO (v=0, 3-2) emission in the disc centered at 344.85 GHz, with 3840 channels, a channel width of 0.42 km s*-1* (effective spectral resolution of 0.82 km s*-1*) and a total bandwidth of 2 GHz.

In all of the scheduling blocks J1337-1257 was used as Bandpass and phase calibrator, with Titan as amplitude and flux calibrator. Calibrations were applied using the pipeline provided by ALMA. The total time on source excluding overheads was 178 min.

3.1 Continuum emission

To study the continuum emission, we use the four spectral windows to reach the highest sensitivity as no CO emission is present in the data (this is discussed below). Figure 2 shows the continuum image using the task CLEAN in CASA 4.4 (McMullin et al., 2007) with natural weights and correcting for the primary beam — Note that the noise increases towards the edges of the image as the primary beam sensitivity decreases. At the center of the image we achieve a rms noise level of 16 $\mu$Jy beam*-1*, which increases to 32 $\mu$Jy beam*-1* at $7\aas@@fstack{\prime\prime}5$. The beam size is $1\aas@@fstack{\prime\prime}1\times 0\aas@@fstack{\prime\prime}7$ with a position angle (PA) of $-70^{\circ}$. In the image three compact sources are detected: 61 Vir’s stellar emission at the center with a total flux of $374\pm 16\ \mu$Jy, which is $2.4\sigma$ higher than the $320\pm 16$ $\mu$Jy predicted photospheric emission assuming a spectral index of -2, thus, it could be due to chromospheric emission at this wavelength (e.g., Loukitcheva et al., 2004; Fontenla et al., 2007); and two other sources to the north of the star with offsets of $4\aas@@fstack{\prime\prime}5$ and $12\aas@@fstack{\prime\prime}5$, and peak fluxes of $360\pm 20\ \mu$Jy and $850\pm 70\ \mu$Jy, respectively. The latter is resolved with a total flux of $2.2\pm 0.3$ mJy within a $2\arcsec$ radius circumference. These two sources are almost certainly the background galaxies previously reported in Wyatt et al. (2012) and their position is overlayed with pentagon markers and labelled as B1 and B2. We also overlay the position of a third background source (B3) detected at 5 GHz with the VLA and not present in the ALMA data. The latter is the southern component of a double-lobed structure with the northern component outside the ALMA primary beam. At 1.4 GHz, B1 was marginally resolved and found to be extended in the north-south direction with a fitted FWHM of $33\arcsec$, therefore, consistent as being the two lobes resolved at 5 GHz.

Although there is no disc emission above $3\sigma$ in the ALMA CLEAN image, significant signal is present in the real component of the visibilities after subtracting the three compact sources. By de-projecting the observed visibilities assuming a disc PA and inclination of $65\degr$ and $77\degr$, respectively (consistent with the Herschel observations, Wyatt et al., 2012), we recover disc emission in the short baselines ($\lesssim 10$ k$\lambda$, see Figure 3), corresponding to extended emission ($\gtrsim 20\arcsec$ or 150 AU). We also overlay the model visibilities of a disc with a flux of 4 mJy and extending from 30 to 140 AU, consistent with the data (see Sec. 4). The imaginary part of the visibilities is consistent with pure noise around zero, which is expected for an axisymmetric centered disc.

We can also recover the disc emission in the image space by integrating the flux inside ellipses of different semi-major axes (with the same PA and aspect ratio or inclination as the disc resolved by Herschel). In this integration we also exclude a $30\degr$ wide wedge in the direction of B1. The resulting radial profile is presented in the top panel of Figure 4. Within $10\arcsec$ the total disc and stellar emission is only $0.8\pm 0.2$ mJy, 2.2 $\sigma$ lower than the derived flux from SCUBA2. If we subtract the stellar emission, the disc is marginally detected at $2.2\sigma$ with a total flux of $0.43\pm 0.2$ mJy. The lower ALMA disc flux could be produced by spatial filtering in the ALMA data due to a lack of short baselines, as the maximum recoverable scale is $6\arcsec$ given the range of baselines in the data. This is illustrated in Figure 3 and demonstrated in Sec. 4, where we fit and simulate the observed visibilities and the SCUBA2 image using a parametric disc model that we use to constrain the disc flux and disc surface density.

We search for any spatially resolved disc emission by azimuthally averaging a CLEAN image of the ALMA data spatially smoothing the emission tapering the visibilities with the Fourier transform of a Gaussian of FWHM of $1\aas@@fstack{\prime\prime}5$. This process degrades the CLEAN beam to a size of $1\aas@@fstack{\prime\prime}8\times 1\aas@@fstack{\prime\prime}5$ increasing the Signal-to-noise ratio (S/N) for extended emission. The azimuthal averaging method also takes into account the disc inclination and PA and is done in wedges of $\pm 30\degr$ along the major axis of the disc. At each radius, the uncertainty is computed based on the uncertainty on each pixel and the number of independent measurements, estimated to be equal to the length of the arc over which we are averaging, divided by the beam’s semi-major axis. The azimuthally averaged intensity is presented in the bottom panel of Figure 4. This shows a marginal disc detection of $0.04\pm 0.01$ mJy beam*-1* at $5\aas@@fstack{\prime\prime}5\pm 0\aas@@fstack{\prime\prime}9$ (where the positional uncertainty is roughly estimated as $\sim$beam semi-major axis/$\sqrt{S/N}$), equivalent to $47\pm 8$ AU, consistent with the inner disc radius constrained to be between $30-40$ AU, depending on the disc model assumed to fit the Herschel observations (Wyatt et al., 2012). Moreover, positive emission, but not significantly above zero apart from the peak at $5\aas@@fstack{\prime\prime}5$, is present from the stellar position to a distance of $11\arcsec$. This is consistent with the positive total flux described before, in other words, with the $2.2\sigma$ detection integrated over all radii.

3.2 CO

Although CO gas of secondary origin has been found in a few young bright debris discs, probably released in collisions of icy solids (e.g., 49 Ceti, $\beta$ Pic, HD 131835, HD 181327 and Fomalhaut, Zuckerman et al., 1995; Dent et al., 2014; Moór et al., 2015b; Marino et al., 2016a; Matrà et al., 2017, submitted) no CO (v=0, J=3-2) emission was detected in 61 Vir ALMA data. Integrating the continuum subtracted channel maps from 30 to 100 AU and radial velocities (RV) in the range $\pm 5.1$ km s*-1with respect to the stellar RV (expected Doppler shift due to Keplerian rotation at 30 AU), we derive an integrated noise level of 27 mJy km s-1*. We can use this to place a $3\sigma$ upper limit to any CO present in the disc. As shown by Matrà et al. (2015) non-local thermodynamic equilibrium (non-LTE) effects can be significant in the low density environments of debris discs; therefore, it is necessary to consider the effect of different gas kinetic temperatures and collisional partner densities — assumed to be electrons released from carbon ionization after the CO gas is photodissociated (as predicted by thermodynamic models, e.g., Kral et al., 2016). Using the tools developed by Matrà et al. (2015) we derive a CO gas mass upper limit of $1.4\times 10^{-6}$ M⊕ using the $3\sigma$ upper limit on the CO flux, the assumed disc extent (30-100 AU) and a carbon ionization fraction of 0.5 and a C/CO abundance of 100 (assumed to be equal to those in $\beta$ Pic, Cataldi et al., 2014; Roberge et al., 2000), which fixes the ratio between electron and CO gas number densities in the disc.

Given the short photodissociation timescale of 120 yr, together with the low dust optical depth, and thus, low collisional rates of solids in the disc, we do not expect to detect CO gas being released in collisions of icy planetesimals in this system. For example, if we assume that planetesimals in the disc have a CO mass fraction of 16%, near the maximum fraction that has been observed in Solar System comets (0.3-16%, Mumma & Charnley, 2011) and similar to other systems with detected exocometary gas (Matra2017Fomalhautinprep), we expect only $\lesssim 10^{-9}$ M⊕ of CO gas in the disc coming out from collisions. Greater amounts of CO gas trapped in ices could come out from icy planetesimals closer in if these are scattered into highly eccentric orbits that can cross the H2O or CO2 snow lines within 10 AU, as suggested by recent ALMA observations of $\eta$ Corvi (Marino et al., 2016a), but this is not detected and no evidence of such scattering has been found so far for 61 Vir.

4 Disc modelling

In order to place better constraints on the total disc flux, disc size, inclination and position angle, we fit a parametric disc model to the SCUBA2 image and ALMA visibility data simultaneously. The model consists of a central star surrounded by a dusty disc and two background point sources (B1 and B2) at the position of the maxima in the ALMA image. The fluxes of the star, B1 and B2 are held fixed at their observed values of 0.37, 0.36 and 0.85 mJy, respectively. Note that the B2 is apparently resolved and could be modelled with an extended component, but this has no effect on the fitted parameters and best fit models.

The dusty disc is assumed to be composed of grains formed by astrosilicates (Draine, 2003), amorphous carbon (Li & Greenberg, 1998), and water ice (Li & Greenberg, 1998), with mass fractions of 70%, 15% and 15%, respectively. We mix the optical constants using the Bruggeman rule (Bohren & Huffman, 1983) and mass-weighted opacities are computed using the Mie theory code of Bohren & Huffman (1983), assuming a Dohnanyi-like size distribution with a power law index of -3.5 (Dohnanyi, 1969), and minimum grain size of 1 $\mu$m, roughly the blow-out size, and a maximum size of 1 cm. We expect larger grains to be present, but we can neglect their thermal emission at this wavelength. The central star is modelled using a stellar template spectrum with a effective temperature of 5500 K111http://www.stsci.edu/hst/observatory/crds/k93models.html (Kurucz, 1979) and a radius of 1.1 $R_{\odot}$ to fit the stellar emission at 0.86 mm. Then, the dust equilibrium temperature at different radii is computed using RADMC-3D222http://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d/ (Dullemond et al., 2016). The disc surface density varies with radius and is parametrized with a power law function as $r^{\alpha}$ from a minimum radius of 30 AU, extending to $R_{\max}$, which is a free parameter as well as $\alpha$ and the total disc flux, $F_{\mathrm{disc}}$. We maintain $R_{\min}$ fixed at 30 AU (best fit value for a model with a sharp inner edge when fitting the Herschel observations and SED, Wyatt et al., 2012). The vertical mass distribution is assumed to be Gaussian with a standard deviation or scale height H that scales linearly with radius as H$=0.1r$. Synthetic images at 0.86 mm are then produced using RADMC-3D with an inclination, $i$, and PA that are also left as free parameters. In total there are 5 free parameters that we vary to fit the observations.

Model visibilities are computed at the same uv points as the ALMA observations (e.g., Marino et al., 2015, 2016a, 2016b). To speed up the simulation of model visibilities, we average the ALMA data with a time and frequency bin of 90 s and 1.88 GHz, respectively. This averaging is small enough both in time and frequency to ensure that the time and frequency smearing are smaller than $0\aas@@fstack{\prime\prime}1$ ($\ll$ synthesised beam). We simulate the SCUBA2 observation by convolving the model image with a two dimensional Gaussian with a FWHM of $13\arcsec$.

To find the best fit we use a Bayesian approach, sampling the parameter space using the python module emcee, which implements Goodman & Weare’s Affine Invariant MCMC Ensemble sampler (Goodman & Weare, 2010; Foreman-Mackey et al., 2013). The posterior distribution is defined as the product between the likelihood function and our prior distributions. The first is defined as $\exp(-\chi^{2}/2)$, with $\chi^{2}=\chi^{2}_{\rm ALMA}+\chi^{2}_{\rm SCUBA2}$, with

[TABLE]

where the sum goes over the uv points of the previously averaged visibilities, $V_{\mathrm{data,i}}$. The estimated error $\delta V_{\mathrm{data,i}}$ is calculated based on the intrinsic dispersion of the visibilities over one scan with the task statwt from casa 4.7. On the other hand, $\chi^{2}_{\rm SCUBA2}$ is defined as the squared sum over every pixel of the difference between the SCUBA2 and model image (convolved with the $13\arcsec$ beam), divided by the pixel rms. The pixel rms is empirically estimated by measuring the dispersion on the unsmoothed SCUBA image that has uncorrelated pixel noise.

The prior probabilities of the parameters are assumed to be uniform. We restrict $R_{\max}$ to be between 30 and 250 AU, $\alpha$ from $-5$ to 5, $F_{\mathrm{disc}}$>0, PA from $0^{\circ}$ to $90^{\circ}$ and $i$ from $45^{\circ}$ to $90^{\circ}$ (priors based on the previous Herschel observations).

To demonstrate that there is disc emission in the ALMA data that can be better constrained by adding the SCUBA2 image to the fitting process, in Figure 5 we present the marginalised distributions of $i$ and PA when fitting only the ALMA data and constraining $R_{\max}$ to values below 140 AU as any disc emission beyond that would lie outside the ALMA primary beam. Even though disc emission above 3$\sigma$ is not present in the reconstructed ALMA image (see Figure 2), but only when integrating the emission, we find that the disc orientation can still be constrained and matches with the previous estimates from Herschel observations (blue lines).

Figure 6 presents the marginalised distributions of $R_{\max}$, $\alpha$ and $F_{\mathrm{disc}}$, when ALMA visibilities and the SCUBA2 image are combined in the analysis. The disc orientation is better constrained, with PA$=59\pm 5$ and $i=82\pm 4^{\circ}$, consistent with the Herschel observations (PA$=65^{\circ}$ and $i=77^{\circ}$), and within the limits obtained from fitting the ALMA data alone (see Figure 5). Regarding the disc structure, we find that $\alpha$ peaks at zero on its marginalised posterior distribution and is constrained between -0.2 and 3.5 (68% confidence), but still consistent within the 95% confidence region with the value of -1 (see Figure 6) found by fitting the Herschel observations, which was also poorly constrained (Wyatt et al., 2012). If we restrict $i$ between $70-80^{\circ}$ (using the prior information from Herschel images), we can improve our constraints on the slope, finding $\alpha=0.1^{+1.1}_{-0.8}$. Therefore, we conclude that the surface density distribution is not very centrally concentrated.

For example, we can discard a scattered disc that has an initial characteristic surface density proportional to $r^{-3.5}$ (e.g., Duncan & Levison, 1997). The collisional evolution of such an scattered disc has been studied analytically by Wyatt et al. (2010). We find that for 61 Vir parameters, i.e. assuming $t=4.6Gyr$, $\alpha=2.5$, pericentre at 30 AU and $M_{\rm disc}\sim 10^{-2}-10^{2}$ M⊕, the resulting surface density should be significantly peaked at 30 AU (pericentre) and decrease steeply with radii, inconsistent with our observations (see their Figure 5). Although the analytic model used by Wyatt et al. (2010) could overestimate the surface density of dust at low radii as it is the case for low eccentricities.

On the other hand, a flat distribution could be expected in the context of an extended disc with a wide range of semi-major axes and small eccentricities, collisionally evolved after being stirred (e.g., Schüppler et al., 2016; Geiler & Krivov, 2017, see Sec. 5). We also find that $R_{\max}$ is peaked at $\sim 150$ AU, consistent with the maximum radius of at least 100 AU derived with Herschel. However, if $\alpha<0.5$ then the maximum radius is not well constrained as the surface brightness decreases with radius ($B(r)\propto r^{\alpha-0.5}$). $F_{\mathrm{disc}}$ peaks above zero ($3.4\sigma$), and is constrained to be $3.7^{+1.2}_{-1.1}$ mJy; however, this is highly dependent on $R_{\max}$ and $\alpha$. For example, if $\alpha\sim-1$ then $F_{\mathrm{disc}}<4$ mJy (95% confidence).

We also try to vary $R_{\min}$ and leave it as a free parameter, but we find that it is not well constrained in these observations. The posterior marginalised distribution of $R_{\min}$ is close to flat with a peak at the inner boundary set to 5 AU. With a smaller $R_{\min}$ the disc surface brightness decreases which fits best the ALMA visibilities, while conserving the total flux to fit the SCUBA observations. Therefore, we decide to leave $R_{\min}$ fixed based on the previous Herschel and SED information that are inconsistent with $R_{\min}\ll 30$ AU.

In Figure 7 we compare simulated observations of different models and their residuals when subtracted from the real observations. The first column shows the best fit model from the posterior distribution presented above with a total flux of 3.5 mJy, $\alpha=0$, $R_{\max}=150$ AU, PA$=65^{\circ}$ and $i=77^{\circ}$, which has a reduced chi-squared $\chi^{2}_{\mathrm{red}}$=1.0028838 ($N\sim 6\times 10^{6}$). The second column shows a model with $\alpha=-1.0$, $R_{\max}=250$ AU and $F_{\mathrm{disc}}=2$ mJy, i.e. the most likely disc flux for this $\alpha$. This model is still consistent with having no disc emission above $3\sigma$ in the reconstructed image and has $\chi^{2}_{\mathrm{red}}$=1.0028840 ($1.4\sigma$ difference with the first model). The third column corresponds to a model similar to the second, but with a less extended disc with $R_{\max}=90$ AU and $F_{\mathrm{disc}}=3$ mJy, increasing the surface brightness of the disc to levels above $3\sigma$ in the simulated observation (Figure 7f), which translates to significant negative residuals (Figure 7i) and $\chi^{2}_{\mathrm{red}}$=1.0028897 ($35\sigma$ difference with the previous model). We also find that the image reconstruction suffers from flux loss due to an insufficient number of short baselines and the size of the primary beam ($17\arcsec$). For the models in the 1st, 2nd and 3rd columns, we recover integrated fluxes of 0.6, 1.0 and $2.3\pm 0.2$ mJy, respectively. From the best fit values of the parameters (i.e. $R_{\max}\gtrsim 150$ AU), the SCUBA2 measured flux ($5.0\pm 1.2$ mJy) and the simulated observations which show that a compact disc would be detectable, we conclude that the disc of planetesimals must be broad and not concentrated in a single or a few narrow rings, which could not have been resolved by Herschel (model 3 in Sec. 1).

5 Steady state collisionally evolved disc model

It is generally assumed in debris discs that the surface density of millimetre-sized grains can be simply scaled to derive the distribution of the total solid mass in discs. This is true under the assumption that the size distribution from big to small bodies remains fixed. However, using detailed numerical simulations with the Analysis of Collisional Evolution (ACE) code (Krivov et al., 2006), Schüppler et al. (2016) recently showed that the radial profile of the vertical optical depth can deviate considerably from the distribution of planetesimals, when considering this more realistic model of how the grain size distribution evolves at different radii. For example, when assuming a maximum planetesimal size of 100 and 200 km in diameter, they found that the optical depth (dominated by the smallest grains in the disc) stays roughly constant as a function of radius between 10 and 100 AU, even though the total surface density decreases with radius. This effect is not due to radiation pressure affecting small dust grains, but due to the evolution of the size distribution at different radii. Specifically, the difference arises when the largest planetesimals in the disc (that dominate the disc mass) are not collisionally evolved, but the smallest grains are already in collisional equilibrium. Using a three phase analytic model for the size distribution, Geiler & Krivov (2017) confirmed this effect and explored how it changes depending on the primordial conditions of the disc.

This implies that even if we assume that the primordial distribution of solids in a debris disc is close to a standard Minimum Mass Solar Nebula (MMSN) with a radial distribution with an exponent of -1.5 after the protoplanetary disc disperses (Weidenschilling, 1977a; Hayashi, 1981), or any model for the initial surface density profile of an accreting protoplanetary disc (e.g., Kuchner, 2004; Raymond et al., 2005; Chiang & Laughlin, 2013), the radial distribution of dust grains with lifetimes shorter than the age of the system could have a significantly different radial dependence. Therefore, the surface density exponent for millimetre grains derived in Sec. 4 cannot be simply extrapolated to the total surface density of solids in 61 Vir.

Here, we aim to study the expected surface density of millimetre grains in a broad debris disc undergoing collisional evolution, and how that depends on the choice of maximum planetesimal size. We do this by using a simple numerical prescription that simulates the size distribution using size bins and assuming that the size distribution is in quasi steady state. This means that the mass loss rate due to catastrophic collisions in each size bin is balanced by the input from fragmentation of larger bodies in destructive collisions, which inputs mass into the bin. The maximum size in collisional equilibrium, $D_{c}$, corresponds to the one having a collisional lifetime equal to the age of the system. This method is described in detail in Wyatt et al. (2011) (see Sec. 2.4.2, 2.5 and 2.9 therein) and can reproduce the morphology (slope and wiggles) seen in more detailed numerical simulations (e.g. using the ACE code, Löhne et al., 2008).

Our model is composed of a 1 M⊙ star at the center and a debris disc spanning 1 to 300 AU. The primordial mass surface density of solids is assumed to be that of a MMSN: $\Sigma_{0}(r)=(r/{\rm 1\ AU})^{-1.5}$ M⊕ AU*-2*, with an initial size distribution of solids proportional to $D^{-3.7}$, though the main results presented below are independent of this choice. The minimum size of solids in the cascade is set to 0.8 $\mu$m, which is the blow-out size assuming a star of 1 L⊙ and 1 M⊙, and an internal density of solids of 2700 kg m*-3*. Grains smaller than this are immediately lost from the disc. We explore different maximum diameters ($D_{\max}$) between 1-100 km. The disc is assumed to be pre-stirred or stirred on a timescale much shorter than the age of the system, i.e. initially having velocities high enough so collisions between planetesimals are destructive and result in a collisional cascade. This is accounted by setting the mean eccentricity ($e$) and inclination ($I$) of the particles to be 0.05 and 1.4*∘* ($e/2$), respectively, which defines the relative velocities of the particles. These velocities are calculated as $v_{\rm rel}=v_{\rm K}(1.25e^{2}+I^{2})^{1/2}$ (valid for Rayleigh distributions of $e$ and $I$, Lissauer, 1993; Wetherill & Stewart, 1993), where $v_{\rm K}$ is the Keplerian velocity on a circular orbit. Hence, the relative or impact velocities are a 6% of $v_{\rm K}$.

Furthermore, in our model destructive collision are only caused by impactors with specific energies greater than the disruption threshold or planetesimal strength ($Q_{\rm D}^{\star}$), which depends both on the size and impact velocity. The disruption threshold has been studied in laboratory experiments (e.g., Fujiwara et al., 1989; Davis & Ryan, 1990; Ryan et al., 1991) and with numerical simulations of colliding basalt and icy bodies (e.g., Benz & Asphaug, 1999). It is well known that for small bodies bound by cohesive binding forces, $Q_{\rm D}^{\star}$ decreases with size up to the size where self-gravity becomes important, and then $Q_{\rm D}^{\star}$ increases with size. Therefore, we assume the following prescription

[TABLE]

where $Q_{\rm D,s}$, $Q_{\rm D,g}$, $b_{\rm s}$ and $b_{\rm g}$ are parameters that depend on the specific composition of solids in the disc. The dependence on the relative or impact velocity is inspired by the results from Stewart & Leinhardt (2009). We use $Q_{\rm D,s}=500$ J kg*-1*, $Q_{\rm D,g}=0.03$ J kg*-1*, $b_{\rm s}=-0.37$, $b_{\rm g}=1.36$ and $v_{0}=3$ km s*-1* values consistent with Basalt in simulations from Benz & Asphaug (1999). The choice of Basalt is not important for the results presented below. Using the values estimated for planetesimals composed of ice from the same study, we obtain similar results. Finally, we assume a “redistribution function” for the fragments created in a destructive collision proportional to $D^{-3.5}$, with the largest fragment having half the mass of the original disrupted body. The specific dependence on $D$ does not change our results presented below.

We divide the disc in different independent annuli, each one with a total mass of $2\pi r\Delta r\Sigma_{0}(r)$, with $\Delta r=2er$, which fixes the initial total mass in each radial bin. At a given radius, we solve for the steady state size distribution by equating the mass loss rate and gain in each size bin that is smaller than $D_{c}$, the largest object that is in collisional equilibrium. The mass in size bins larger than $D_{c}$ is held fixed to the primordial distribution as they have lifetimes or collisional timescales longer than the age of the system and have not had enough time to significantly evolve. The timescale to reach quasi steady state or damp perturbations is the same as the collisional timescale; therefore, our quasi steady state assumption is valid for sizes smaller than $D_{c}$. To find the specific $D_{c}$ we solve for the steady state size distribution varying $D_{c}$, until finding the specific size bin with a lifetime equal to the age of the system (or with a difference smaller than a 10%). In the resulting size distribution bins for planetesimals larger than $D_{c}$ retain their original masses, while the masses in all smaller bins are anchored to $D_{c}$ and their size distribution set by the collisional equilibrium condition. As the system age increases, $D_{c}$ increases, and the size distribution evolves, and thus, the total and mm-sized dust mass too.

If $D_{c}>D_{\max}$, i.e. the lifetime of the biggest planetesimal is shorter than the age of the system ($t_{\rm age}$), the mass in every bin is scaled as

[TABLE]

where $t_{\rm age}$ is the age of the system, and $M^{\prime}_{0}$ is the mass distribution in collisional equilibrium when $D_{c}=D_{\max}$, or when the system had an age equal to the lifetime of the biggest planetesimal, $t_{\mathrm{c}}(0)$. Equation 3 is valid if the mass loss rate is proportional to M2, which is the case in our models as the collisional lifetime is inversely proportional to the mass in the cascade. The evolution of the surface density of solids at 10 AU is illustrated in the top panel of Figure 8. The main relevant feature of this evolution is that when $D_{c}<D_{\max}$, the mass in the small size bins decreases more slowly than it would when the entire size distribution is in equilibrium ($D_{c}=D_{\max}$).

In Figure 9 we present the evolution of three discs varying $D_{\max}$ from 1 to 100 km (top and middle), and changing the stellar mass and luminosity together with $\Sigma_{0}$ and $D_{\max}$ to fit 61 Vir disc properties (bottom panel), i.e. its surface density of mass in mm-sized grains and disc inner edge (see Sec. 5.1). The surface density of the total mass in solids ($\Sigma$, left column) evolves with time similarly to analytic models (e.g., Wyatt et al., 2007), increasing with distance as expected up to the characteristic radius, $r_{c}$, at which the largest planetesimal in the disc has a lifetime equal to the age of the system, i.e. $t_{\mathrm{c}}(0)=t_{\rm age}$, and from there decreasing with radius as $\Sigma_{0}(r)$. This radius depends on the initial total solid mass and on $D_{\max}$ as the three panels in the first column show, with $r_{c}$ being smaller for larger $D_{\max}$ or lower initial mass as the rate of collisions is reduced. The surface density of mass in mm-sized grains ($\Sigma_{\mathrm{mm}}$, right column) behaves in a way similar to the optical depth described in Schüppler et al. (2016) and fractional luminosity in Geiler & Krivov (2017), mimicking $\Sigma(r)$ for $r<r_{c}$, but considerably flatter compared to $\Sigma(r)$ at $r>r_{c}$. This is because $D_{c}<D_{\max}$ and $D_{c}$ decreases with $r$ outside $r_{c}$, so $\Sigma_{\rm mm}$ is less depleted for larger $r$. The net effect is that $\Sigma_{\rm mm}(r)$ is almost constant, even though $\Sigma(r)$ decreases with $r$.

To illustrate the differences in the evolved size distribution at different radii, in the bottom panel of Figure 8 we compare the size distribution at 1 Gyr with $D_{\max}=100$ km and $r=10$, 40 and 100 AU. At 10 AU (red line) the disc evolves fast as relative velocities are higher and all the size bins are in collisional equilibrium. At 40 AU (purple line) relative velocities are slower, hence $D_{c}\sim 20$ km and only smaller bodies are in collisional equilibrium; therefore, the mass in small bodies is highly depleted compared to the primordial, while the total mass in solids has not decreased significantly. At 100 AU (blue line) relative velocities are even slower, $D_{c}\sim 2$ km and the mass in small bodies is less depleted compared to the primordial than at 40 AU. Even though the surface density of the total mass in solids at $t=0$ and 1 Gyr is higher at 40 AU than at 100 AU, the mass surface density in solids smaller than 1 km is approximately the same at both radii after 1 Gyr. This causes the slope of $\Sigma_{\mathrm{mm}}$ to flatten out and be almost constant at large radii where $t_{\mathrm{c}}(0)>t_{\rm age}$ as mentioned above. We also observe a very similar evolution for the vertical optical depth in the disc, consistent with Schüppler et al. (2016).

This behaviour that makes $\Sigma_{\rm mm}$ to be almost flat can be understood analytically if we consider a planetesimal strength approximated by two broken power laws and a continuous size distribution with three regimes: i) small bodies in collisional equilibrium with a size distribution proportional to $D^{-q_{1}}$; ii) large bodies with gravity dominated strengths in collisional equilibrium with a size distribution proportional to $D^{-q_{2}}$; and iii) largest planetesimals with lifetimes longer than the age of the system that conserve their primordial size distribution proportional to $D^{-q_{3}}$. The value of $q_{1}$ and $q_{2}$ are strictly related to the dependence on $D$ for $Q_{\rm D}^{\star}$, with (Durda & Dermott, 1997; O’Brien & Greenberg, 2003)

[TABLE]

where $b_{i}$ is the slope or exponent of $Q_{d}^{\star}$ in the strength or gravity-dominated regime. Therefore, assuming reasonable values for $b_{s}$, $b_{g}$ and $q_{3}$ we can find an analytic expression for the size distribution at different radii (e.g., Löhne et al., 2008). Moreover, assuming an initial surface density or mass distribution in the disc, we can derive an expression for the fractional luminosity as a function of radius, as shown by Shannon & Wu (2011) (Equation A10 therein) and Geiler & Krivov (2017) (Equation A11 therein). As the fractional luminosity is proportional to the Surface density of small grains, we can rewrite Equation A10 in Shannon & Wu (2011) to find

[TABLE]

where $\Sigma_{0}(r)$ is the primordial total surface density of solids, and $k_{2}$ is equal to $\frac{6-q_{2}}{q_{2}-1}$ and represents the size scaling of the minimum impactor size to cause a catastrophic collision. The expression above is only valid when $D_{c}$ is less than $D_{\max}$, but large enough so it is in the gravity dominated regime ($D_{c}\gtrsim 100$ m). Assuming $q_{1}$=3.6, $q_{2}$=3.0, $k_{2}=1.5$ (values consistent with $b_{s}$ and $b_{g}$ used above), $q_{3}$=3.7 and $\Sigma_{0}(r)=\Sigma_{0}(r/1\ {\rm AU})^{\alpha}$ we find

[TABLE]

Therefore, for $\alpha=-1.5$, $\Sigma_{\rm mm}$ is independent of radius, which matches with our more detailed numerical simulation. Moreover, the flatter $\Sigma_{\rm mm}$ in the evolved size distribution compared to the primordial distribution is independent of $\alpha$ as Equation 6 shows; although a steeper primordial surface density of solids decreasing with radius would result in a steeper surface density of millimetre sized grains with a slope of $0.6\alpha+0.9$. For $q_{3}=3.5$ and 3.9 we still find a flat slope for $\Sigma_{\rm mm}$ of -0.3 and 0.2, respectively. From the results in our simulations we can estimate the dependence of $\Sigma_{\rm mm}$ on $t,D_{\max}$ and $\Sigma_{0}$, by assuming a power law dependence and fitting it to our numerical results. Coupling these with the dependence on $r$ from Equation 6 (only valid for $q_{1}$=3.6, $q_{2}$=3.0 and $q_{3}$=3.7), we find

[TABLE]

where $\Sigma_{0}$ is the initial surface density of solids at 1 AU in units of the MMSN. The factor 2 and the exponents of $-0.4$, $-0.1$ and $0.6$ are the results from a fit to the numerical simulations. Equation 7 is only valid for $r>r_{c}$ and $D_{c}\gtrsim 100$ m. Using Equation A5 in Shannon & Wu (2011) we can also estimate how $r_{c}$ varies with time and $\Sigma_{0}$. Moreover, from our simulations we can derive a dependence on $D_{\max}$ fitting a power law. We find

[TABLE]

Assuming a specific dependence of planetesimal strength on size, equations A10 from Shannon & Wu (2011), and 7 and 8 from this work, together with the disc model presented above, can be used to retrieve the primordial radial distribution of solids from ALMA observations of extended discs if the biggest planetesimals in the disc still conserve their primordial size distribution. Moreover, they can be used to constrain the initial total mass in the disc and maximum planetesimal size. So far there is no evidence of extended debris discs at millimetre wavelengths with a steep slope decreasing with radius (or non consistent with being flat, e.g., Booth et al., 2016); however, even with ALMA (the most sensitive instrument at millimetre wavelengths) the detection and study of extended debris discs is only possible around a few of the brightest systems.

It is worth noting that the maximum planetesimal size in a disc could vary with radius by orders of magnitude as the growth timescales for planetesimals are a steep function of radius and the surface density in solids (e.g., Kenyon & Bromley, 2008). Moreover, stirring could have stopped the growth at different epochs for different radii. Although in our models we assume that the maximum planetesimal size is independent of radius, our prediction for $\Sigma_{\rm mm}(r>r_{c})$ in Equation 7 is not very sensitive to $D_{\max}$. Therefore, our predictions are reasonably valid even if the maximum planetesimal size decreases with radius (as expected in planet formation models). This is already illustrated in Figure 9. If we consider a disc with $D_{\max}$ decreasing from 100 to 1 km between 40-300 AU, then the resulting $\Sigma_{\rm mm}(r)$ at 1 Gyr would be almost the same as the red line in the middle right panel on that Figure, because $\Sigma_{\rm mm}(300\ {\rm AU})$ increases only by a factor of 2 when decreasing $D_{\max}$ from 100 to 1 km. This is due to two opposite effects: 1) for a constant total mass in solids, reducing $D_{\max}$ increases the mass in millimetre sized dust; and 2) reducing $D_{\max}$ makes the collisional evolution faster which reduces the mass in every bin in collisional equilibrium. A similar effect would be present at $r<r_{c}$ making the surface density slope flatter. The maximum planetesimal size is only significantly important to determine $r_{c}$. The opposite scenario, and less likely, in which $D_{\max}$ increases with radius would result in a slightly steeper slope for both $r<r_{c}$ and $r>r_{c}$.

Other differences relative to our assumptions could also change the slope of the millimetre surface density, such as the epoch of stirring (in our simulations we consider a pre-stirred disc), or the mean eccentricity and inclination of particles in the disc, or even the disruption threshold of planetesimals and dust if their composition varies with radius. For example, a different $Q_{\rm D}^{\star}$ would modify the size distribution, changing the slope of the predicted millimetre surface density as Equation 5 shows.

5.1 Application to 61 Vir

In Sec. 4 we find that the observations in the millimetre are best fitted with a disc extending to $\sim 150$ AU, an integrated flux of $3.7\pm 1.2$ mJy and a flat surface density distribution, equivalent to a dust mass of $\sim 2\times 10^{-8}$ M⊕ AU*-2*. In addition, the minimum radius derived from a best fit model of a collisionally evolved disc to the Herschel observations is $\sim 40$ AU. Using the same model for the collisional evolution of a disc described above (replacing the stellar mass and luminosity with 0.88 M⊙ and 0.84 L⊙, Sousa et al., 2008; Wyatt et al., 2012) we find a best match with a primordial surface density between 20-100 times less dense than the MMSN and a maximum planetesimal size between $5-20$ km. These two parameters determine that $t_{c}(0)=4.6$ Gyr at $\sim 40$ AU and $\Sigma_{\mathrm{mm}}(r)\sim 2\times 10^{-8}$ M⊕ AU*-2* for $r>40$ AU.

The need for a low primordial surface density and a maximum planetesimal size of 10 km is due to the low mass in millimetre grains, which scales roughly as $D_{\max}^{-0.1}\Sigma_{0}^{0.6}$ (see Equation 7), together with a large $r_{c}$, that scales roughly as $D_{\max}^{-0.5}\Sigma_{0}^{0.33}$ (see Equation 8). Therefore, we need a very low $\Sigma_{0}$ to fit the millimetre surface density and a low $D_{\max}$ to have $r_{c}\sim 40$ AU given the low $\Sigma_{0}$. From the size distribution we can also determine a vertical optical depth of $2\times 10^{-4}$, a few times higher than the optical depth from Herschel observations and SED fitting, but still consistent considering all the assumptions made in both modelling efforts. For example, a more detailed treatment of radiation pressure could change the value of the optical depth by a factor of a few. The derived maximum planetesimal size and primordial surface density go in the same direction as the ones from Wyatt et al. (2012): the primordial surface density of solids in the disc was much lower compared to the MMSN and with a maximum planetesimal size not much larger than 10 km.

6 Discussion

6.1 A depleted broad disc of planetesimals

In Sec. 4 we found that the debris disc in 61 Vir is broad, extending from 30 to 150 AU or larger radii. If the emission were concentrated in a few $\lesssim$20 AU wide rings of planetesimals the disc would have been detected above 3 sigma in the ALMA map. Moreover, the 2.2$\sigma$ difference between the flux measured by SCUBA2 and ALMA is indicative that there is flux loss in the reconstructed ALMA image due to the disc emission being mostly in structures on scales larger than $6\arcsec$ (50 AU); and thus not recoverable by the range of baselines in the ALMA data. This was corroborated using simulated observations of different broad disc models. Therefore, we conclude that the planetesimal disc must be extended with a wide range of semi-major axes. A different scenario with a population of highly eccentric planetesimals with a small range of semi-major axes is discarded as the derived surface density is flatter than expected in a scattered disc scenario (e.g., Duncan & Levison, 1997) and while collisional erosion can flatten this distribution by preferentially eroding the inner regions this cannot completely erase the density enhancement at the inner edge of the disk (Wyatt et al., 2010).

The inner edge of the disc could be defined by the collisional evolution that has been ongoing for 4.6 Gyr as assumed in Sec. 5, or instead the disc could have been truncated by a yet unseen planet. In the first scenario, the observed inner edge of the disc (30-40 AU) can be explained by a maximum planetesimal size of only about 10 km and primordial surface density of solids 50 times lower than a MMSN. One explanation for why the planetesimals did not grow to larger sizes could be the low surface density of solids which slows down the growth timescales (Kenyon & Bromley, 2008), but could also be because the planetesimals were stirred by a planet closer in hindering their growth.

In the second scenario, in which the inner edge of the disc has been truncated by a planet, the maximum planetesimal size is no longer constrained to be of the order of $\sim 10$ km. However, even if we consider a maximum planetesimal size of 1000 km, the mass of the primordial disc still needs to be a factor of $\sim 10$ lower compared to a MMSN in order to fit the flat surface density of millimetre grains derived in this paper and the Herschel observations (Wyatt et al., 2012). This depletion could arise from the protoplanetary disc phase if the disc had a low mass, or a low efficiency of planetesimal formation, or due to radial drift of solid particles during that gas rich phase which concentrated most of the solid mass in the inner regions (Whipple, 1973; Weidenschilling, 1977b). The radial drift of solids could have also contributed to an in situ formation of the 2-3 planets found within 1 AU of the Star (e.g., Hansen & Murray, 2012).

A variant on the second scenario involves the 30 AU truncation radius being caused by a planet which is no longer present. For example, if the close-in planets formed further out (just inside 30au) and then migrated to their current location accreting and scattering planetesimals on their way in, the early evolution of these close-in planets could be responsible for both the truncation of the outer disc and its stirring (e.g., Alibert et al., 2006; Terquem & Papaloizou, 2007; Kennedy & Kenyon, 2008; Payne et al., 2009; Ida & Lin, 2010).

6.2 Stirring by a yet unseen planet

If 61 Vir b and c formed in situ, then something else must have stirred the disc as these are too far in and not massive enough to have stirred the disc at large radii within 4.6 Gyr (Wyatt et al., 2012). Hence, we propose that an unseen planet at a larger distance and within the 30 AU disc inner edge stirred the disc. Similar to Moór et al. (2015a), using Eq. 6 from Mustill & Wyatt (2009) (valid for planets with eccentricities $\lesssim 0.3$) we can derive lower limits on the eccentricity of such a planet depending on its semi-major axis and mass so the timescale of stirring is shorter than the age of the system. Moreover, the eccentricity imposed on the planetesimals ($e_{f}$) must be higher than a certain value so that their relatives velocities are high enough to cause destructive collisions ($v_{\mathrm{rel,max}}\sim 2e_{\mathrm{f}}v_{\mathrm{K}}$). Here we impose that the forced eccentricity (Equation 8 in Mustill & Wyatt, 2009) must be higher than 0.01 so planetesimals of 5 km diameter undergo destructive collisions with planetesimals of the same size at 150 AU. This is illustrated in Figure 10. The minimum eccentricity decreases with increasing semi-major axis and planet mass as the timescale for stirring is held fixed at 4.6 Gyr. The forced eccentricity must be $>0.01$, which results in a kink in the 0.1 contour (because $e_{f}$ is independent of mass). All other contours are set by the stirring time set equal to the age of the system.

We can add additional constrains if we require planets with a pericenter that does not get closer than 5 mutual Hill radii (see Eq. 9 in Pearce & Wyatt, 2014) from the apocentre of 61 Vir c (a=0.22 AU, e=0.14), i.e

[TABLE]

where $R_{\rm H,q}$ is the Hill radii at pericentre and $a_{\rm c}$ and $e_{\rm c}$ are the semi-major axis and eccentricity of 61 Vir c. In addition, the apocentre of the hypothetical planet d has to be such that it does not get closer than 5 Hill radii to the disc inner edge at $\sim 30$ AU. These two constraints exclude the grey shaded area. Because lower mass planets have higher eccentricity, the maximum semi-major axis decreases with decreasing planet mass for M${}_{\rm plt}\lesssim 10$ M⊕, but also decreases with increasing planet mass as the $R_{\rm H,q}$ gets larger. Finally, using upper limits from RV data from HARPS we can exclude planets more massive than the red line (Wyatt et al., 2012; Kennedy et al., 2015).

With these limits on Mplt and $a_{\rm plt}$ we can conclude that if an unseen planet interior to the debris disc is responsible for stirring the planetesimals up to 150 AU, and has an eccentricity lower than 0.1, then it must be more massive than 10 M⊕ and have a semi-major axis between 10-20 AU. Less massive planets and closer in ($a_{\rm plt}=4-20$ AU) could have stirred the disc, but with $e\gg 0.1$. For the allowed combinations of $M_{\mathrm{plt}}$ and $a_{\mathrm{plt}}$ even a highly eccentric planet will not induce an eccentricity higher than the observed on 61 Vir b and c or cause close encounters (see Figure 5 in Read & Wyatt, 2016). Moreover, an eccentric planet will impose an eccentricity on the disk which may be detectable by imaging of the disc (Wyatt et al., 1999). While there is no evidence of any asymmetry, the constraints are weak, both because the imaging is in the far-IR where the transition from pericentre to apocentre glow occurs (Pan et al., 2016), and because the disc would look symmetric if the pericentre is aligned with the minor axis of the disc projected in the sky.

The equations used to derive the minimum eccentricity are only valid for $e_{\rm plt}\lesssim 0.3$. Planet eccentricities higher than 0.3 could be overestimated as the predicted stirring timescales are longer than expected for $e>0.3$ (see Figure 1 in Mustill & Wyatt, 2009). Therefore, the lower limits presented in Figure 10 are only representative of the constraints expected for $e_{\rm plt}\lesssim 0.3$.

6.3 Background sources

As noted before by Wyatt et al. (2012), none of the detected compact sources (B1 and B2) are co-moving with 61 Vir, therefore, we can assume these are background objects. B1 together with B3 and its northern counterpart (which lies outside the ALMA primary beam) are probably related to an active galactic nucleus (AGN, Condon et al., 1998). B1 at the center is consistent with compact emission from dust heated by an AGN, or with flat-spectrum synchrotron emission typical of a radio galaxy core, or with dust associated with a nuclear starburst, or some combination of these possibilities. On the other hand, the two lobes are consistent with synchrotron emission that we do not expect to detect in the sub-millimetre (sub-mm), given their steep spectra.

The most northern source detected by ALMA (B2) was previously detected by Herschel from 160 to 500 $\mu$m and is resolved by the ALMA synthetic beam with a size of $\sim 2\arcsec$. This is larger than expected for the z$\sim$1-3 sub-mm galaxy population (Smail et al., 1997, e.g.,), where ALMA has measured typical sub-mm sizes of $\ll 1\arcsec$ (e.g., Ikarashi et al., 2015; Simpson et al., 2015), slightly smaller than their typical radio sizes (Biggs et al., 2011), but is consistent with a dusty starburst at a rather lower redshift, as we would suspect from its relatively blue PACS/SPIRE SED.

7 Summary and Conclusions

We have presented the first resolved millimetre study of 61 Vir, a planetary system with two confirmed RV planets within 1 AU and a debris disc at tens of AU. Combining ALMA and SCUBA2/JCMT observations we found that at 0.86 mm the total disc emission is $3.7\pm 1.2$ mJy, the disc extends from 30 to at least 150 AU, and has surface density exponent of millimetre grains of $0.1^{+1.1}_{-0.8}$. This implies that the parent planetesimal disc is broad with a wide range of semi-major axes. The alternative scenario of a highly scattered disc with planetesimals with a common pericentre is discarded given the constraints on the surface brightness of the disc. No CO gas emission was detected in the disc, although even if planetesimals are rich in CO and releasing gas through collisions, we predict that any emission should be below our detection limit.

We developed a full disc collisional evolution model based on previous numerical work that can reproduce some of the results obtained in more detailed simulations, but in a much more computationally efficient approach. These models predict that the surface density of millimetre grains and optical depth radial profiles do not necessary match with the surface density of the parent bodies, tending to be flatter in regions of the disc where the age of the system is shorter than the collisional lifetime of the biggest planetesimals. This can be used to constrain the primordial surface density distribution of solids and maximum planetesimal size for extended discs for reasonable assumptions on the eccentricity, inclination and strength of planetesimals. For example, with this model we can reproduce the observations if 61 Vir debris disc started with a surface density $\sim 50$ times more depleted in solids compared to a MMSN, and with planetesimals that did not grow more than 5-20 km in size so the disc is collisionally depleted at $r<40$ AU. However, these conclusions are based on the assumption that the inner edge of the observed disc is set by the collisional evolution of the disc. If instead the inner edge is set by other mechanism, e.g. planet-disc interaction, then the maximum planetesimal size is no longer constrained, but the primordial surface density would still need to be depleted by a factor of $\sim 10$ compared to the MMSN.

Finally we discussed and constrained the mass, semi-major axis and eccentricity of a planet stirring the disc located between the known RV planets and the inner edge of the disc. We found that in order to have stirred the disc out to 150 AU, the planet must be more massive than 10 $M_{\oplus}$ and a semi-major axis between 10 and 20 AU if it has an eccentricity lower than 0.1. Otherwise, for higher eccentricities it could have a lower mass and a semi-major axis between 4 and 20 AU.

Acknowledgements

We thank Pablo Roman and Simon Casassus for providing us the tool uvsim to simulate model visibilities. We also thank Matthew J. Read for useful discussion. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00359.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and 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. MCW, LM and AS acknowledge the support of the European Union through ERC grant number 279973. GMK is supported by the Royal Society as a Royal Society University Research Fellow. AS is partially supported by funding from the Center for Exoplanets and Habitable Worlds. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium.