ALMA resolves CI emission from the β Pictoris debris disk
Gianni Cataldi, Alexis Brandeker, Yanqin Wu, Christine Chen, William Dent, Bernard L. de Vries, Inga Kamp, René Liseau, Göran Olofsson, Eric Pantin, Aki Roberge
DDraft version July 2, 2018
Typeset using L A TEX default style in AASTeX61
ALMA RESOLVES C I EMISSION FROM THE β PICTORIS DEBRIS DISK
Gianni Cataldi,
1, 2, 3, 4, ∗ Alexis Brandeker,
3, 4
Yanqin Wu, Christine Chen,
6, 7
William Dent, Bernard L. de Vries,
3, 4, 9
Inga Kamp, Ren´e Liseau, G¨oran Olofsson,
3, 4
Eric Pantin, and Aki Roberge Subaru Telescope, National Astronomical Observatory of Japan, 650 North Aohoku Place, Hilo, HI 96720, USA Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, Konkoly-Thege Mikl´os ´ut15–17, 1121 Budapest, Hungary AlbaNova University Centre, Stockholm University, Department of Astronomy, Stockholm, Sweden Stockholm University Astrobiology Centre, Stockholm, Sweden Department of Astronomy & Astrophysics, University of Toronto, Ontario, Canada Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA Department of Physics and Astronomy, The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD, 21218, USA ALMA/ESO, Alonso de Cordova 3107, Santiago, Chile Scientific Support Office, Directorate of Science, European Space Research and Technology Centre (ESA/ESTEC), Keplerlaan 1, 2201 AZNoordwijk, The Netherlands Kapteyn Astronomical Institute, University of Groningen, The Netherlands Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden Laboratoire AIM, UMR 7158, CEA/DSM—CNRS—Universit´e Paris Diderot, IRFU/SAp, 91191 Gif sur Yvette, France Exoplanets and Stellar Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
ABSTRACTThe debris disk around β Pictoris is known to contain gas. Previous ALMA observations revealed a CO belt at ∼
85 au with a distinct clump, interpreted as a location of enhanced gas production. Photodissociation converts COinto C and O within ∼
50 years. We resolve C I emission at 492 GHz using ALMA and study its spatial distribution.C I shows the same clump as seen for CO. This is surprising, as C is expected to quickly spread in azimuth. We derivea low C mass (between 5 × − and 3 . × − M ⊕ ), indicating that gas production started only recently (within ∼ (cid:38) M Moon . Keywords: circumstellar matter — stars: individual ( β Pictoris) — submillimeter: planetary systems— radiative transfer — techniques: interferometric — methods: observational
Corresponding author: Gianni [email protected] ∗ International Research Fellow of Japan Society for the Promotion of Science (Postdoctoral Fellowships for Research in Japan (Standard)). a r X i v : . [ a s t r o - ph . E P ] J un INTRODUCTIONIn debris disk systems, the continuous collisional destruction of larger bodies such as comets or asteroids producesabundant amounts of dust, with the smallest grains quickly removed by radiation pressure (Backman & Paresce1993; Wyatt 2008). A debris disk provides evidence that planetesimal-sized bodies were formed during the earlierprotoplanetary phase (Artymowicz 1997; Matthews et al. 2014) and gives us the opportunity to study the propertiesof the building blocks of planets. Studying these systems is therefore intimately linked to our efforts of understandinghow planets form.The debris disk around the young (23 ± β Pictoris hasbeen used as a laboratory to study the early evolution of planetary systems ever since its discovery by the InfraredAstronomical Satellite (IRAS) (Aumann 1985). Smith & Terrile (1984) obtained the first resolved image showing anedge-on disk. Since then, the properties of the dust disk have been extensively studied with observations at multiplewavelengths. Today, we know that the belt of parent bodies is located at ∼
100 au (Dent et al. 2014) and that β Pichosts a giant planet (e.g. Chilcote et al. 2017, and references therein) with a semi-major axis of ∼
10 au (Lecavelier desEtangs & Vidal-Madjar 2016; Wang et al. 2016), first imaged by Lagrange et al. (2009, 2010).Even before the dust disk was discovered, evidence from optical and ultraviolet (UV) absorption lines suggested thepresence of gas around β Pic (Slettebak 1975; Slettebak & Carpenter 1983). The β Pic disk is thus part of a smallsub-sample of debris disks where gas has been detected. Currently, there are about 20 such gaseous debris disks known(e.g. Redfield 2007; Mo´or et al. 2011; Roberge et al. 2014; Lieman-Sifry et al. 2016; Matr`a et al. 2017a).The spatial distribution of the gas around β Pic has been studied with resolved observations in the optical andrecently with the Atacama Large Millimeter/submillimeter Array (ALMA). These data showed that the gas disk isradially extended (with some species traced out to several hundred au) and in Keplerian rotation (Olofsson et al. 2001;Brandeker et al. 2004; Nilsson et al. 2012; Dent et al. 2014). Besides this stable component, time-variable absorptionfeatures shifted with respect to the systemic velocity are attributed to exocomets evaporating in vicinity of the star(e.g. Vidal-Madjar et al. 1994; Kiefer et al. 2014, and references therein). This latter phenomenon has also been seenaround a number of other (mostly young) A-type stars (e.g. Welsh & Montgomery 2015, and references therein).Similarly to the dust, the gas in the β Pic disk is thought to be continuously produced from the destruction ofsolid material rather than leftover from the protoplanetary phase. Evidence for such a secondary scenario comes forexample from theoretical arguments on the dynamical lifetime of the gas (Fern´andez et al. 2006). Also, models of theexcitation of the CO 3–2 and 2–1 transitions observed by ALMA imply that not enough H is present in the disk toshield CO from photodissociation over the lifetime of the disk, thus the necessity of a gas replenishment mechanism(Matr`a et al. 2017b). Studying this secondary gas opens up the interesting possibility to constrain the composition ofthe parent bodies (e.g. Kral et al. 2016; Matr`a et al. 2017b, 2018).To date, multiple metallic species such as C, O, Na, Al or Ca have been detected (e.g. Brandeker et al. 2004; Robergeet al. 2006; Brandeker et al. 2016). Recently, the first detection of hydrogen was reported by Wilson et al. (2017). COremains the only molecule detected so-far (e.g. Dent et al. 2014; Matr`a et al. 2018). The observed elemental abundancesare strikingly different from solar abundances. While the abundance of H is much lower than solar (Wilson et al. 2017),C is highly overabundant with respect to other metals (Roberge et al. 2006; Cataldi et al. 2014). Fern´andez et al.(2006) showed that the carbon overabundance provides a braking mechanism, preventing other metals that are stronglyaffected by radiation pressure (such as Na) from being quickly ejected from the system. Carbon also plays a crucial rolein determining the excitation conditions of atomic fine-structure or molecular rotational transitions (e.g. Zagorovskyet al. 2010). This is because in a secondary, hydrogen-depleted disk, collisional excitation is likely dominated byelectrons, and ionisation of carbon is the main electron source.Using spectrally resolved observations of C II with Herschel /HIFI, the spatial distribution of carbon was constrainedby Cataldi et al. (2014). They found that most of the carbon gas is located at ∼
100 au, with tentatively more emissionfrom the south-west (SW) side of the disk. At about the same time, ALMA spatially resolved CO 3–2 emission,revealing a clump of emission on the SW side of the disk at ∼
85 au (Dent et al. 2014). The CO clump coincides witha radial peak of the millimetre continuum (Dent et al. 2014) and a clump seen at mid-IR wavelengths (Telesco et al.2005; Li et al. 2012). Dissociation by interstellar UV photons limits the lifetime of CO in the disk to significantly lessthan an orbit (Visser et al. 2009). It is thus clear that CO needs to be produced continuously, and is a source for Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and withimportant participation from NASA. both C and O. This might provide a natural explanation for the observed super-solar abundance of C with respectto metals such as Na (Xie et al. 2013). The CO itself is believed to be derived from the destruction of volatile-rich,cometary bodies, where the clump corresponds to a location of increased collision rate and thus gas production. Severalhypothesis have been put forward to explain the existence of the clump. Firstly, it could be the location of a giantcollision. The clump then results from the fact that the orbits of the collision debris always all go through the collisionpoint (Jackson et al. 2014). Secondly, the clump could be due to cometary bodies trapped in a resonance with anoutward migrating, yet unseen giant planet (Wyatt 2003, 2006). Using ALMA follow-up observations of the CO 2–1transition at higher resolution, Matr`a et al. (2017b) dismissed the giant collision scenario based on the large radialextent of the clump. Thirdly, Nesvold & Kuchner (2015) proposed that the interaction between a spiral density waveand a vertical displacement wave, both induced by β Pic b, can produce an azimuthaly asymmetric collision rate.Assuming that indeed all C and O is derived from the dissociation of CO, Kral et al. (2016) modelled the hydro-dynamical evolution of carbon and oxygen in the disk. They assumed that the produced atomic gas is subject tothe magneto-rotational instability (MRI) and predict an atomic accretion and decretion disk (Xie et al. 2013; Kral &Latter 2016).C I was previously seen in absorption against β Pic (Vidal-Madjar et al. 1994; Jolly et al. 1998; Roberge et al. 2000).Recently, Higuchi et al. (2017) presented the first observation of C I in emission. In this paper, we present the firstspatially resolved observations of C I , revealing its distribution in the disk. Our paper is organised as follows. Wedescribe the observations in section 2 and present the results in section 3. In section 4, we present simple gas emissionmodels to study the total mass and spatial distribution of the carbon gas. We discuss our results in section 5 andconclude in section 6. OBSERVATIONS AND DATA REDUCTIONWe observed the β Pic disk using band 8 receivers of ALMA on December 19, 2015 during the ALMA cycle 3 EarlyScience campaign (project ID 2013.1.00459.S, PI: Brandeker). The observations were split into two execution blocks(EBs). The total integration time was 2.1 h (with 1.2 h on β Pic) and the median precipitable water vapor (pwv) was0.6 mm with standard deviation of 0.1 mm. The array consisted of 36 antennas arranged in a hybrid configurationcontaining both short and long baselines ranging from 15 m to 6.3 km. In principle, we are thus sensitive to angularscales between ∼ (cid:48)(cid:48) and ∼ (cid:48)(cid:48) on the sky . However, the visibilities from the longer baselines were affected by largeatmospheric phase fluctuations and therefore flagged during the calibration process (see below). The observations wereexecuted as a mosaic to obtain uniform sensitivity over the whole disk. One pointing was centred on the star and twoadditional pointings were centred at ± (cid:48)(cid:48) along the position angle of the disk. The size of the primary beam is 11.8 (cid:48)(cid:48) .Antenna elevations varied between 27 ◦ and 60 ◦ .We placed three spectral windows, each with 1920 channels and a channel spacing of 488 kHz (total bandwidth937.5 MHz), onto the following transitions: C I P – P at 492.16 GHz, CS 10–9 at 489.75 GHz and SiO 11–10 at477.50 GHz. For C I , this corresponds to a channel spacing of 0.30 km s − and an effective spectral resolution of0.34 km s − (spectral averaging factor N = 2). In addition, a fourth spectral window with 128 channels and a totalbandwidth of 1875 MHz was placed at 479 GHz in order to observe the dust continuum.The data were calibrated using CASA uvcontsub taskwithin
CASA . We then imaged the visibilities using the
CLEAN task within
CASA . In order to increase our surfacebrightness sensitivity, we applied a taper of 1 (cid:48)(cid:48) , thus significantly reducing the contribution of the remaining long see ALMA Cycle 6 Technical Handbook, section 3.6 (Spatial Filtering), https://almascience.nao.ac.jp/documents-and-tools/cycle6/alma-technical-handbook see ALMA Cycle 6 Technical Handbook, Table 5.2
200 150 100 50 0 50 100 150 200x [AU]0.50.00.51.0 f r a c t i o n o f p e a k
10 5 0 5 10x [arcsec]
Figure 1.
Left : Moment 0 of the C I
492 GHz emission from the β Pic disk. Contours are drawn at intervals of 3 σ , with negativecontours drawn as dotted lines. The beam is illustrated as white ellipse in the lower left. The dashed rectangle illustrates theregion used to measure the total flux. Right : Emission profile along the x-axis of the moment 0 map, normalised to the peakvalue. The gray shaded area illustrates the ± σ error interval. The symbol in the upper left shows the projection of the beamonto the x-axis. baselines (at 492 GHz, an angular scale of 1 (cid:48)(cid:48) corresponds to a baseline length of 126 m). We also produced a continuumimage using CLEAN with the same taper, combining all spectral windows except the one centred onto CS 10–9, whichis in a region of bad atmospheric transmission and therefore particularly noisy.Because of the low antenna elevations and the suboptimal array configuration, the sensitivity of the data is signifi-cantly worse than requested. Consequently, the data initially did not pass Quality Assurance 2 (QA2). However, westill decided to publish the data as the signal-to-noise ratio (SNR) is sufficient to provide new information on the β Picsystem. RESULTS3.1.
Line emission
The CS 10–9 and SiO 11-10 lines remained undetected. We detect and resolve C I P – P emission. Figure 1(left) shows the moment 0 map, produced by integrating the data cube along the spectral axis within ± − (withrespect to the systemic velocity of the star, assumed to be v heliocentric = 20 . − , Brandeker 2011). The figure hasbeen rotated to align the horizontal direction with the main dust disk (position angle of +29 . ◦ , Lagrange et al. 2012).The beam size is 1.18 (cid:48)(cid:48) × (cid:48)(cid:48) (23 au ×
19 au) with a major axis position angle with respect to the North of 44 ◦ . Weachieve a 1 σ sensitivity of ∼ . × − W m − Hz − sr − ( ∼
70 mJy beam − ) in the individual channels. We performphotometry by considering a rectangular aperture extending ±
115 au in the horizontal direction (measured from thestellar position, see Figure 1) and ±
30 au in the vertical direction (measured from the mid-plane). This yields a totalflux of (1 . ± . × − W m − (9 . ± . − ), where the error is random without any systematic calibrationuncertainty taken into account. We did not correct for the primary beam, as it changes the total flux only within thequoted error bars (the same applies for the continuum, section 3.3). The measured flux is consistent with the valueof (1 . ± . × − W m − (10 . ± . − ) derived from single dish observations by Higuchi et al. (2017).The same method yields 3 σ upper limits of 5 . × − W m − (3.6 Jy km s − ) for CS 10–9 and 2 . × − W m − (1.8 Jy km s − ) for SiO 11-10.Figure 1 (right) shows the emission profile along the x-axis of the moment-0 map, obtained by integrating within ±
30 au in the vertical direction. Both the moment-0 map and the emission profile are suggestive of an asymmetrywith a peak on the SW side of the disk. The same asymmetry is seen for the CO emission (Dent et al. 2014; Matr`aet al. 2017b) and tentatively also for C II (Cataldi et al. 2014). However, the SW/NE flux ratio within the rectangularbox of Figure 1 is not significant at 0 . ± . ab by propagating the error like this: σ a/b = b σ a + a b σ b ). The SW/NE ratio of the peaks in the emission profile is 1 . ± .
5. Thus, there is not necessarilymore flux on the SW side of the disk, but the flux seems to be more compact.Figure 2 shows the position-velocity (pv) diagram (i.e. the data cube integrated along the vertical spatial directionwithin ±
30 au). This figure thus shows the radial velocity of the emission as a function of the projected position along
Figure 2.
Position-velocity diagram of the C I emission. Contours are drawn at 3 σ intervals. The spectro-spatial resolutionis illustrated in the lower left by the white rectangle. The two solid lines show the radial velocity (for circular Keplerian orbitsaround a 1.75 M (cid:12) star seen edge-on) at 50 au and 220 au, the approximate radial extent of the CO (Matr`a et al. 2017b). Thedashed lines show the region in x-v-space used to derive a moment 0 map and emission profile with improved SNR. the disk mid-plane. By using the pv diagram, we can constrain the radial distribution of the emission despite theedge-on orientation of the disk. Indeed, figure 2 also shows the radial velocity for circular Keplerian orbits at 50 auand 220 au around a 1.75 M (cid:12) star (Crifo et al. 1997), seen edge-on. This is the approximate radial extent of the CO asfound by Matr`a et al. (2017b). As can be seen, no significant C I emission is detected beyond these lines, suggestingthat most C I emission is confined to the same region as the CO.We may use the pv diagram to define a region in p-v space that contains all significant emission. Then, whenintegrating over the spectral axis, we only take data points within this region into account (rather than everythingwithin ± − as was done to produce Figure 1). The region is illustrated by the dashed lines in Figure 2 andthe resulting moment 0 map and emission profile are shown in Figure 3. The SW/NE flux ratio (within the samebox as in figure 1) is now 1 . ± .
2. Most importantly, the SNR of the emission profile is significantly improved andthe SW/NE asymmetry more clearly visible. We measure a SW/NE peak ratio of 1 . ± .
4. Thus, the significanceof the peak ratio asymmetry is only marginal. However, in Figure 3 (right) we also show the profiles of the CO 3–2and 2–1 emission (see Figure 2 of Matr`a et al. 2017b), for which the SW/NE flux ratios are 1 . ± .
08 (CO 2–1) and1 . ± .
13 (CO 3–2) and the peak ratios are 1 . ± .
14 (CO 2–1) and 1 . ± .
14 (CO 3–2). The C I emission followsthe CO 2–1 emission surprisingly well, with the same distinct peak on the SW side of the disk and the same asymmetricextent (out to ∼
150 au in the NE and ∼
100 au in the SW). This suggest that the asymmetry is indeed real. This is asurprising result. The asymmetry in CO can be readily understood from its short (less than one orbit) lifetime due tophotodissociation: if CO is primarily produced in a clump, photodissociation prevents CO from spreading azimuthally,thus preserving the asymmetry. On the other hand, the C produced from CO photodissociation is expected to spreadin azimuth within a few orbits. In section 5.5 we discuss possible solutions to this puzzle.
200 150 100 50 0 50 100 150 200x [AU]0.50.00.51.0 f r a c t i o n o f p e a k CICO 2-1CO 3-210 5 0 5 10x [arcsec]
Figure 3.
Same as Figure 1, but using only the region in p-v space illustrated by the dashed lines in Figure 2 to improve theSNR. For the moment 0 (left), no contours are shown since by construction, the noise levels are now non-uniform and dependon x . The color scale is the same as in Figure 1. In the emission profile (right) we also included the CO 3–2 and 2–1 emissionas observed by Matr`a et al. (2017b). From the moment 0 maps, it is also apparent that the observed emission is slightly tilted, with the NE side belowand the SW side above the midplane of the main outer disk (defined by z = 0). A similar tilt is seen for CO (Dentet al. 2014; Matr`a et al. 2017b). As is discussed by Matr`a et al. (2017b), two reasons might be imagined for thisobservation: either the gas disk is indeed tilted with respect to the dust disk, or the tilt is a projection effect, resultingfrom an azimuthally asymmetric gas disk in combination with a slight deviation from a perfect edge-on inclination.Interestingly, the inner dust disk seen in scattered light has a similar tilt (e.g. Milli et al. 2014; Apai et al. 2015;Millar-Blanchaer et al. 2015). This dust component known as warp or secondary disk is localised inwards of 80 AU(Lagrange et al. 2012) and seems thus slightly inwards of the gas. Also, it has already been suggested that this innerdust disk is not perfectly edge-on (e.g. Milli et al. 2014; Millar-Blanchaer et al. 2015).The reader interested in the procedures to estimate the errors quoted in this section is referred to Appendix A.3.2. Deprojection of the C I emission Assuming that the gas follows circular Keplerian orbits, we can use the pv diagram (figure 2) to obtain a face-onview of the emission (e.g. Dent et al. 2014; Matr`a et al. 2017b). Indeed, each point in the pv diagram corresponds totwo points in the xy-plane of the disk (where we define y as the coordinate running along the line of sight), one infront and one behind the sky plane. There remains the degeneracy of how to distribute the flux of a given pv pointamong the two points in the xy-plane. As discussed by Matr`a et al. (2017b), the degeneracy can be broken if the diskis not perfectly edge-on. However, for simplicity, we follow Dent et al. (2014) and assume that the disk is edge-on. Ourprimary interest is to illustrate the radial extent of the emission and the position of the clump. We thus choose to placeall flux in front of the sky plane, but other physically motivated choices are possible (Dent et al. 2014). Figure 4 showsthe deprojection. Points of the pv diagram with | x | <
40 au were masked because the radial velocity in this regionstends towards zero for all radii, i.e. it becomes difficult to assign a radius to a given radial velocity. Emission is seenapproximately out to 150 au. The clump appears at a similar position angle as seen in CO (Dent et al. 2014). Detailsof the deprojection procedure are described in appendix C. Since the optical depth of the emission is not negligible,the deprojection does not show the true distribution of the emission in the xy -plane but rather how much emissionthe observer receives from various locations in the xy -plane.3.3. Continuum
The continuum image at ∼
485 GHz is shown in figure 5. The 1 σ noise level is 4 . × − W m − Hz − sr − (1.3 mJy beam − ). The beam size is 1.19 (cid:48)(cid:48) × (cid:48)(cid:48) (23 au ×
19 au) with the position angle of the major axis being46 ◦ . We measure a total flux of (1 . ± . × − W m − Hz − (112 ± ±
140 aufrom the star along x and ±
30 au from the mid-plane along z ) indicated in the figure (see appendix A for details onthe error calculation). The measured flux is consistent with a Rayleigh-Jeans extrapolation of the flux measured by Figure 4.
Deprojection (face-on view) of the C I emission derived from the pv diagram (figure 2), where we chose to place allflux in front of the sky plane (the line of sight runs in the positive direction of the y -axis). Points with | x | <
40 au are masked.
ALMA at 870 µ m (Dent et al. 2014), and is a factor of ∼ . ± .
14. In anycase, the asymmetry is analogous to what is seen in thermal mid-IR images between 8.7 and 18.3 µ m (Telesco et al.2005) as well as for C and CO. MODELLING4.1.
Simple estimation of the total carbon mass
In this section, we estimate the total carbon mass in the β Pic disk under some simplifying assumptions. With a 1Dmodel, we want to reproduce the C I flux measured in the present work and the C II flux measured with Herschel /PACS(Pilbratt et al. 2010; Poglitsch et al. 2010) by Brandeker et al. (2016, OBSID 1342198171) (we sum the flux values ofthe PACS spaxels listed in their table 1).
Herschel /HIFI also measured the C II flux (Cataldi et al. 2014). However,the flux measured by PACS is likely a better estimate of the total C II flux from the disk, because the HIFI half powerbeam width is only ∼
200 au.We need to compute the statistical equilibrium of the level populations, as in general the emission cannot be assumedto be in local thermal equilibrium (LTE). The gas in the β Pic disk is thought to be of secondary origin and thusdepleted in hydrogen (e.g. Matr`a et al. 2017b). We thus assume that the dominant collider species is electrons, andthat photoionisation of carbon is the main electron source (Fern´andez et al. 2006; Kral et al. 2016), i.e. the density ofionised carbon equals the electron density. Under these assumptions, the total carbon mass is estimated in two steps.For a given kinetic temperature, we first determine the amount of ionised carbon necessary to reproduce the C II Figure 5.
ALMA continuum image of the β Pic disk at 480 GHz. Contours are drawn at intervals of 3 σ . The white dashedrectangle indicates the region within which the flux was measured. emission. The second step is to determine the amount of neutral carbon necessary to reproduce the C I flux observedby ALMA, using the electron density derived in the first step.To solve the statistical equilibrium and radiative transfer, we use pythonradex , a python implementation of the RADEX code (van der Tak et al. 2007) with atomic data from the LAMDA database (Sch¨oier et al. 2005). This codeuses an escape probability formalism to solve the radiative transfer and assumes the geometry of a uniform sphere.In general, the radiative transfer depends on the geometry of the emitting region. The gas observed around β Pic isclearly non-spherical and non-symmetric. The assumption of a uniform sphere is thus a clear limitation of this model.The masses derived here should therefore be considered first order estimates.We choose the size of the sphere such that its volume corresponds to the volume of an elliptic torus with semi-majoraxis of 35 au in the radial direction and semi-minor axis of 10 au in the vertical direction (see Figures 1,2 and 3). Notethat even in the optically thin case, such an assumption on the scale over which the emission is produced is necessary,unless the emission is in LTE. This is because for a given mass, the electron density depends on the assumed volume.We include radiative excitation and de-excitation (hereafter simply (de-)excitation) by the CMB, the stellar radiation(at 85 au) and the dust continuum, where the latter is taken at 85 au as seen in figure B1 of Kral et al. (2017) and isthe dominant component. However, it turns out that including radiative (de-)excitation from these sources does notchange our results because (de-)excitation is dominated by collisions and spontaneous emission.We then consider a wide range of kinetic temperatures T kin from 40 K to 1000 K. For lower temperatures, the C II line quickly becomes strongly optically thick such that meaningful constraints on the mass are no longer possible (i.e.the model cannot reproduce the observed flux regardless of how much the mass is increased). However, based on ourmore detailed modelling in section 4.2, we deem lower temperatures unlikely. Detailed thermodynamical modellingby Kral et al. (2016) also suggests that the temperature is above ∼
50 K within 100 au, although in their model, thetemperature drops to 20 K at 150 au.Figure 6 shows the derived total carbon mass as a function of T kin . The figure also shows the individual C andC + masses. To assess the importance of non-LTE effects, we also show corresponding curves for which LTE has beenassumed. Both C I and C II are generally speaking close to LTE. For higher temperatures, the C mass required toreproduce the observed C I flux is higher in LTE than in non-LTE. This simply happens because in LTE, higher levelsget populated more quickly with increasing temperature, thus de-populating the level that produces the C I emission.The maximum optical depths encountered for the considered temperature range are 0.5 for C I and 3.9 for C II . https://github.com/gica3618/pythonradex http://home.strw.leidenuniv.nl/~moldata/ From Figure 6, we determine a total carbon mass between 5 × − and 3 . × − M ⊕ . The lower bound is quiterobust to changes of the size of the emitting region as it is close to the LTE value. For example, increasing or decreasingthe radius of the emitting sphere by 50% does not change the lower bound by more than 15%. On the other hand, itis clear that the upper bound is more uncertain as it can quickly increase if one allows for lower temperatures and/orsmaller emitting volumes (i.e. increased optical depth). Another parameter that can strongly affect the optical depth(and thus the upper bound of the mass range) is the assumed line broadening parameter. Here we used b = 0 .
57 km s − ,derived in appendix B.Previous studies estimating the carbon gas mass in the β Pic disk had only the C II flux and line profile atdisposal. These more detailed models derived higher masses by using the spectrally resolved C II observations fromHerschel/HIFI: Cataldi et al. (2014) obtained 1 . × − M ⊕ while Kral et al. (2016) derived 1 . × − M ⊕ (the latterstudy also used an upper limit on the C I flux). In the Kral et al. (2016) model, the total C mass is dominated byionised carbon that is located beyond ∼
100 au. The C I flux is of little importance for the total C mass budget ofthat model. The discrepancy can thus not be explained by the fact that we include the C I measurement. On theother hand, most of the carbon in the Cataldi et al. (2014) model is located between 150 and 300 au with an ionisationfraction of roughly 50%, i.e. there is a significant contribution of neutral carbon to the total mass. However, ourALMA data show no C I emission beyond ∼
120 au. That model indeed over-predicted the C I flux by a factor of ∼ II line profile (see their figure 2b) clearly suggests that C II emission beyond 150 au is needed to fitthe line core (this is also an issue for our 3D models (section 5.2) that generally do a bad job in fitting the C II linecore). So maybe there is largely ionised carbon gas beyond ∼
150 au present and the reason why we derive a lowerionised carbon mass here is because we assumed (based on the C I data) a too small volume (i.e. too high density andthus more excitation, and therefore less mass needed). However, even when increasing the volume of our 1D model,we still derive a lower mass of ionised carbon compared to Kral et al. (2016) and Cataldi et al. (2014) and can thusnot fully resolve the discrepancy, which might be due to the different assumptions of the models. The relatively smallHIFI beam (FWHM of ∼
200 au) could also play a role. For example, deviations from axisymmetry in the distributionof ionised carbon could introduce additional uncertainty in the mass estimates of Cataldi et al. (2014) and Kral et al.(2016) since C II emission from large radii is only detected by HIFI if it arises close to the line of sight towards thestar. 4.2.
3D modelling
In this section, we model the 3D spatial distribution of the carbon gas using our new ALMA observations of C I and the previously published spectrally resolved Herschel /HIFI observation of C II (Cataldi et al. 2014, OBSID1342238190). For a given, arbitrary distribution of carbon gas, we first solve the ionisation balance in each grid cell.We use the ionisation rate from Zagorovsky et al. (2010) (scaled with distance to the star) and assume that the gas isoptically thin to ionising photons. The star is the main source of ionising photons, but ionisation from the interstellarradiation field (ISRF) is also included. The ionisation balance is solved analytically by assuming that all electronscome from the photoionisation of carbon. The ionisation fraction thus only depends on the distance to the star, thelocal gas density and the kinetic temperature (via the recombination coefficient). Next, we use the derived electrondensity to solve the statistical equilibrium (SE) of the level populations using atomic data from the LAMDA database.The following processes can (de)excite an atom: spontaneous emission, collisions with electrons (we neglect othercolliders) and radiative (de)excitation by line photons and the background radiation field, where the latter is assumedto be composed of the CMB, the star and the dust continuum. For the dust continuum, we employ the field shown infigure B1 of Kral et al. (2017) (that figure actually shows the field in the midplane of the disk at the position angleof the clump (L. Matr`a, private communication), but for simplicity, we only take the radial dependence into account).In principle, a full radiative transfer computation is necessary to solve the SE. However, we take a simplified approachand assume that the line emission is essentially optically thin. Basically, while the emission can become opticallythick along the disk mid-plane (in our best-fit models, the maximum optical depth is typically of the order of 1), it isoptically thin in the vertical direction, i.e. the photon escape fraction is high. Thus, we assume that the backgroundfield is not attenuated by the gas (all atoms are subject to the full background field) and that (de)excitation by linephotons can be neglected. These assumptions make the SE easy to solve. We further discuss this approximation inappendix D. The background radiation turns out to be unimportant for our models.0 T kin [K]10 -3 m a ss [ M ⊕ ] total CC C + Figure 6.
Mass of C , C + and total C (derived from the model described in section 4.1) as a function of the assumed kinetictemperature. For the full lines, the statistical equilibrium has been solved, while for the dashed lines, LTE has been assumed. Having solved the SE, we compute the emitted spectrum for each grid cell, red- or blue-shifting it according to itsradial velocity. The final step is then to ray-trace the emission along the line of sight to take optical depth into account.The result is a model data cube that can be compared to observations. For simplicity, we consider isothermal modelsand fix T kin = 75 K everywhere. We found that for lower temperatures, the C I /C II flux ratio tends to be too high,while the inverse is true for higher temperatures.To compare to the C II data, we take the corresponding model data cube, multiply by the HIFI beam and integratespatially to get a model HIFI spectrum (Cataldi et al. 2014). For the C I ALMA observations, we convolve the modeldata cube to the same spatial and spectral resolution as the observations and multiply by the primary beam. Theresidual moment 0 maps (respectively pv diagrams) shown in figures 7, 8 and 11 are obtained by subtracting the modelmoment 0 map (pv diagram) from the data moment 0 map (pv diagram) shown in figure 1 (figure 2).4.2.1.
Uniform ring model
We first consider a simple, symmetric model consisting of a ring with uniform surface density. The number densityreads n ring ( r, z ) = Σ √ πH ( r ) · exp (cid:16) − z H ( r )) (cid:17) if r min ≤ r ≤ r max r and z are cylindrical coordinates (with z perpendicular to the disk mid-plane), H ( r ) is the scale height, r min and r max define the radial extent of the ring and Σ is the constant surface density. The scale height in hydrostaticequilibrium for a vertically isothermal disk is given by (e.g. Armitage 2009) H ( r ) = (cid:115) kT r µm p GM ∗ (2)1 Table 1.
Explored parameter space and best fit values for the ‘ring’ and ‘ring + clump’ models fitting the C I ALMA dataand the C II Herschel /HIFI data simultaneously. The number of values explored for each parameter is denoted by n . We alsoindicate whether the values are uniformly distributed in linear or logarithmic space. M ring and M clump are the masses of thering and the clump respectively. The reference mass of β Pic is assumed M β Pic = 1 . M ∗ (Crifo et al. 1997). We also listthe constant surface density of the ring Σ, the mid-plane number density of the ring n mid at 85 au and, for the ‘ring + clump’model, the number density at the centre of the clump n clump resulting from the superposition of the ring and the clump.parameter unit min max n spacing best fit‘ring’ ‘ring + clump’ M ∗ M β Pic M ring M ⊕ . × − . × − . × − . × − r min au 30 70 3 lin 50 50 r max au 120 160 3 lin 120 120 M clump M ⊕ . × − . × − . × − r c au 70 100 3 lin 70 σ xy au 20 40 3 lin 40Σ cm − . × . × n mid ( r = 85au) cm −
140 140 n clump cm − with k the Boltzmann constant, T the temperature, m p the proton mass, G the gravitational constant and M ∗ thestellar mass. For the mean molecular weigh µ , we follow Matr`a et al. (2017b) and assume µ = 14. At 85 au and for T = 75 K, the scale height is 4.2 au.We compute a grid of models over the parameter space listed in the first four rows of table 1. Note that we also testdifferent values for the (dynamical) stellar mass that affects the orbital velocities (and scaleheight). To each model,we assign a χ measure by using expressions analogous to equation 2 in Booth et al. (2017) (we take the correlationof neighbouring data points into account by using an appropriate noise correlation ratio for each data set). We sumthe χ from the Herschel /HIFI data (C II ) and the ALMA data (C I ).Figure 7 shows the C I residuals in the xz and xv plane as well as the C II line emission of the model with thelowest χ . The corresponding model parameters are given in table 1. The model provides a decent fit to the data, butunsurprisingly is not able to reproduce the clump observed in the SW. Thus, we consider a more complicated modelby adding a clump to the uniform ring in the next section. Such a model has no direct physical justification, but isuseful to empirically constrain the spatial distribution of the gas and get an estimate of the gas mass.4.2.2. Uniform ring with a clump
The clump is modelled as n clump ( x, y, z ) = n · exp (cid:18) − ( x − x ) + ( y − y ) σ (cid:19) · exp (cid:18) − z H ( r )) (cid:19) (3)where x runs along the disk mid-plane in the sky plane and y along the line of sight (and x + y = r ). Furthermore, x = r c cos( φ c ) and y = r c sin( φ c ) designate the center of the clump, where we have defined r c as the clump’s radialdistance to the star and φ c as its azimuthal angle in the x-y-plane. We fix φ c = − ◦ (Matr`a et al. 2017b). Thestandard deviation of the clump density distribution in the x - y -plane is σ xy and the density at the centre n . Thetotal carbon number density is then given by n C = n ring + n clump . We compute models over the parameter space listedin table 1. Figure 8 is the analogue of figure 7 for the best ‘ring + clump’ model. As can be seen, the model does aslightly better job in modelling the clump. 4.2.3. Eccentric gas distribution
The ‘ring + clump’ model of the previous section is purely empirical and does not have a direct physical justification.As mentioned earlier, proposed mechanisms to get such a morphology are either a giant collision or resonant trapping2
Figure 7.
Comparison of the uniform ring model to the data. The top row shows the residual moment 0 map (left) and pvdiagram (right) for the C I emission. Contours are in steps of 3 σ . The bottom row compares the modeled C II emission (redlines) to the HIFI data (black lines). of cometary bodies by a migrating giant planet. However, as we discuss in section 5.5, based on our new C I data, wedeem both of these possibilities unlikely.Another way to get a morphology with an emission clump on one side of the disk is to relax the assumption ofcircular orbits and instead consider gas on eccentric orbits. In section 5.5.6, we discuss how such orbits could arise inthe first place. As an example, we here consider a model where the eccentricity of the orbits is uniformly distributedbetween two values e min and e max and where all orbits share the same pericentre and the same argument of periapsis.The pericentre is then a region of higher density and can mimic a clump as seen in the observations.The derivation of the gas density for this model is given in appendix E. We again compute a grid of models withthe parameters listed in table 2. Figure 9 shows a face-on view of the mid-plane carbon gas density and figure 10 thepv diagram of the modelled C I emission. As is seen in figure 11, this model is similarly effective in reproducing theobserved C I emission, although it has some difficulties to reproduce the C II line. We emphasize that the purpose ofthis model is merely to demonstrate that eccentric gas orbits are able to reproduce the general morphology with theclump. A more detailed model will be presented in a forthcoming paper.An interesting consequence of eccentric orbits is that the gas along the line of sight towards the star has non-zero radial velocity. Thus, emission (or absorption) towards the star is shifted with respect to the systemic velocity.For example, our best-fit model predicts that the emission peak towards the star appears 0.4 km s − blue-shifted.Another consequence is an additional velocity broadening compared to the circular case, with broadening parameter b ecc ≈ − (this is smaller than and thus consistent with the broadening measured from the pv diagram, seeAppendix B). The velocity shift and broadening depend on the eccentricity and orientation of the disk. To verify thisvelocity shift and constrain the model, we would need to know the absolute stellar radial velocity to better accuracythan ∼
80 m s − (for a 5 σ -detection of the shift).3 Figure 8.
Same as figure 7, but for the ‘ring + clump’ model.
Table 2.
Same as table 1, but for the parameter space explored by the models with eccentric orbits. The total mass of themodel is M tot , the pericentre distance is r per and the argument of pericentre is ω .parameter unit min max n spacing best fit M ∗ M β Pic M tot M ⊕ . × − . × − . × − e min - 0 0.4 5 lin 0 e max - 0.1 0.5 5 lin 0.3 r per au 60 100 3 lin 80 ω deg -135 135 7 lin -45 The column densities of neutral carbon towards the star for the three models discussed in sections 4.2.1, 4.2.2 and4.2.3 are (6–7) × cm − , slightly above the value of (2–4) × cm − measured in absorption by Roberge et al.(2000). For ionised carbon, the models range between 5 × cm − and 1 . × cm − , while Roberge et al. (2006)report 2 × cm − . 4.3. Absence of an atomic accretion disk
Recently, Kral et al. (2016) presented a model that computes the temperature, ionisation and hydrodynamicalevolution of the atomic carbon and oxygen in the β Pic disk. In this model, the atomic gas is produced in a parentbelt from photodissociation of CO and then evolves viscously under the influence of the MRI, i.e. it forms an accretionand decretion disk. Kral et al. (2016) predict that the carbon is mostly neutral in the inner parts of the accretion disk.4
Figure 9.
Face-on view (mid-plane number density) of the C gas for the best fit model with eccentric orbits described in section4.2.3.
However, figures 2 and 4 indicate that no atomic accretion disk has formed (yet) as there seems to be little C I emission inside of ∼
50 au (but see also section 5.4). To confirm this, we compute the C I emission expected from themodel and compare it to our data. We thus take the Kral et al. (2016) distribution of neutral carbon and electrons,temperature profile and scale height (see their figure 9) and compute the C I emission with the methods described insection 4.2. Figure 12 shows the moment 0 map and pv-diagram of the model and the residuals when subtracting themodel from the ALMA C I data. It is clear from this figure that the accretion profile predicted by Kral et al. (2016)produces too much C I emission in the inner regions of the disk. This is also clearly visible in the pv diagram wheretoo much emission is present at high velocities. DISCUSSION5.1.
Dynamical mass of the star
For all the 3D models described in section 4.2, the grid searches indicate that a dynamical mass of 0 . M ∗ (with M ∗ = 1 . M (cid:12) the assumed stellar mass, Crifo et al. 1997) is preferred to reproduce the data. While this result hasto be interpreted with care given that we did not derive any error bars, it is interesting to note that Olofsson et al.(2001) derived the same 0 . M ∗ dynamical mass by modelling spatially and spectrally resolved Na I emission. Theyascribed their finding to radiation pressure opposing the gravity of the star. In fact, radiation pressure on Na is sostrong that a braking mechanism is required to keep it on the observed Keplerian orbits and preventing it from beingblown out (e.g. Liseau 2003; Brandeker et al. 2004). Fern´andez et al. (2006) suggested that all species affected byradiation pressure are largely ionised and coupled into a single fluid by Coulomb interactions, with carbon acting as abraking agent. In this situation, neutrals are also expected to be coupled (the only exception would be neutrals thatdo not ionise yet still feel a significant radiation force). Thus, one might actually expect that Na I and C I share thesame dynamics. The dynamical mass (respectively the radiation pressure coefficient) is coupled to the composition ofthe gas (Fern´andez et al. 2006) and could thus give interesting information on the elemental abundances. However,5 Figure 10.
PV diagram of the C I emission for the best fit model with eccentric orbits, degraded to the spectro-spatialresolution of the ALMA data. given the unknown uncertainty of the derived dynamical mass which could also be influenced by the assumptions ofour models, we do not draw any further conclusions at this stage.5.2. Issue with the C II line profile core All the 3D models presented in section 4.2 have difficulties to reproduce the core of the C II line profile (see figures7, 8 and 11). The issue is particularly pronounced for the model with eccentric orbits. The fits to the C II line profileby Cataldi et al. (2014) indicate that the line core requires ionised carbon beyond ∼
150 au to be present. As alreadydiscussed in section 4.1, there might thus exist a gas component with a high ionisation fraction beyond ∼
150 au. Thiscomponent might also be required to prevent the metals observed there from being blown out by radiation pressure (seesection 5.7). It is possible that our simple models are unable to capture this extended gas component. For example,a more complex density and/or temperature profile might be required to reproduce it. Note that while the Cataldiet al. (2014) model fits the C II line profile, it would be a bad fit to our new C I data as it strongly overpredicts theC I flux. 5.3. Overall timescale of CO and C production
The strong spatial correlation between the C I and CO emission (Figure 3) suggests a scenario where the carbonis mainly produced from photodissociation of CO, i.e. the mass loss rate of CO equals the production rate of C.Thus, from the CO mass loss rate and our estimate of the total C mass, we can estimate the time since C (and CO)production started in the β Pic disk, provided that no carbon has yet been removed.5.3.1.
Revised CO lifetime
It was previously thought that photodissociation of CO in the β Pic disk is dominated by the ISRF. For example,taking self-shielding into account, Matr`a et al. (2017b) calculated a CO lifetime of ∼
300 a in the clump against theISRF. Here we revise this value, showing that photodissociation from the star actually dominates over the ISRF.6
Figure 11.
Same as figure 7, but for the best fit model with eccentric orbits.
We compute the CO photodissociation rate using photodissociation cross sections from the Leiden Observatorydatabase of ‘photodissociation and photoionisation of astrophysically relevant molecules’ (Visser et al. 2009; Heayset al. 2017). We calculate an unshielded lifetime against the ISRF (Draine 1978; Lee 1984) of ∼
130 a.For the stellar spectrum, the basis is a PHOENIX model as described in Fern´andez et al. (2006). This model is thencomplemented with data in the UV from the Space Telescope Imaging Spectrograph (STIS) aboard the
Hubble SpaceTelescope (HST) (Roberge et al. 2000) as well as from the the Far Ultraviolet Spectroscopic Explorer (FUSE) (Bouretet al. 2002; Roberge et al. 2006). This is important because β Pic shows additional emission in the UV above thepredictions from a standard stellar atmosphere model. These additional UV photons impact the calculated lifetimeof CO. The lifetime corresponding to the observed stellar spectrum is ∼
70 a at 85 au. However, since the light weobserve has travelled through the full CO and C column, this is actually the lifetime of a CO molecule sitting behindthis column. To obtain the unshielded lifetime against the star, we have to multiply by the CO self-shielding function(table 6 in Visser et al. 2009) and the shielding function of the C ionisation continuum (Rollins & Rawlings 2012),evaluated at the full CO and C column densities against the star (Roberge et al. 2000). Thus, the unshielded lifetimeagainst stellar photons at 85 au is ∼
20 a.From Matr`a et al. (2017b) and Roberge et al. (2000), we know the vertical column density at the clump locationand the horizontal column density of CO against the star respectively. By applying a rough scaling based on thedeprojected CO emission in Matr`a et al. (2017b), we estimate horizontal and vertical column densities at the clumplocation and along the line of sight to the star. For C , the horizontal column density against the star is taken fromRoberge et al. (2000), while for the vertical column density, we take our best-fit ‘ring + clump’ model as a reference,which is also used to scale the Roberge et al. (2000) horizontal column density to the clump location. For the shieldingtowards the star, it is not difficult to see that the average CO molecule experiences half of the total column density http://home.strw.leidenuniv.nl/~ewine/photo/ Figure 12.
Model C I emission (left column) computed from the predictions by Kral et al. (2016), and residuals (right column)when subtracting the model from the ALMA data. The top row show the moment 0, and the bottom row the pv diagram.Contours are in steps of 3 σ , with negative contours drawn as dotted lines. against the star. For a Gaussian sphere, the average column density per molecule towards the ISRF turns out to bearound half of the column density seen from the centre of the sphere. We choose shielding factors corresponding tothese average column densities. Applying these shielding factors to the unshielded lifetimes derived above, we findthat the overall CO lifetime at 85 au is similar in the clump and the disk: ∼
50 a.5.3.2.
C production timescale
Using the CO mass 3 . × − M ⊕ derived by Matr`a et al. (2017b), a CO lifetime of 50 a leads to a CO mass lossrate of 1 . × kg s − . Under the assumptions that 1) the CO mass loss rate is constant, 2) C is produced only fromCO photodissociation (C might also be produced from the destruction of other carbon-bearing molecules such as forexample methane, but we expect this to be a minor contribution given typical abundances in Solar System comets(e.g. Mumma & Charnley 2011)) and 3) no C is removed from the system, we can use the carbon mass derived insection 4.1 to calculate the time-scale over which CO and C production has been ongoing: t CO = t C = N C / ˙ N CO with N C the total number of carbon atoms and ˙ N CO the destruction rate of CO in molecules per second. Using the C massrange derived from our 1D model (section 4.1), we find a timescale between 1 . × and 1 . × years. When usingthe mass from the best fitting ‘ring + clump’ model of section 4.2.2, we find a timescale of ∼ density. Removing carbon by radiation pressure seems unlikely. First, bothneutral and ionised C do not feel any radiation pressure around β Pic. However, as shown by Fern´andez et al. (2006),ions are coupled into a single fluid via Coulomb interactions. But the effective radiation pressure on this fluid is notsufficient for blow out (this explains why e.g. Na is seen in Keplerian rotation despite feeling strong radiation pressure),so C is not blown out as part of the fluid either. Recondensation onto dust grains is also irrelevant (Grigorieva et al.2007, although they considered the recondensation of water, but similar arguments apply for C gas). Finally, accretionby a planet seems unlikely. Gap-opening planets (Jupiter-class) are the best candidates. But even such planets havea finite accretion efficiency (typically 75%–90%) limited by the leakage of flow across their gaps (Lubow & D’Angelo2006). So to increase the estimated lifetime by a factor of 10 or more, the hypothetical planet would have to sit closeto the observed belt and accrete with an efficiency of 90% or higher. In addition, planets down to 1 M
Jup have beenexcluded beyond 30 au by direct imaging searches (Absil et al. 2013).5.4.
Why is there no accretion disk?
Kral et al. (2016) predicted the formation of a C accretion disk based on thermo-hydrodynamical modelling. However,we showed in section 4.3 that no such accretion disk is present. If the C gas production indeed started as recently asestimated in section 5.3, the absence of an accretion disk is actually not surprising. Indeed, the timescale for viscousaccretion from a radius r can be estimated by t vis = r αc s H (4)with c s the sound speed and H the scale height. We assume c s = (cid:113) kTµm p with T = 75 K and µ = 14, and H as inequation 2 with r = 85 au. A very high α (cid:38) α , the gas has not yet had enough time to form an accretion disk.However, a certain amount of carbon is still needed in the inner regions of the disk where metals such as Na or Feare seen in Keplerian rotation. These species are strongly affected by radiation pressure and C is needed to preventthem from being blown out. As was shown by Fern´andez et al. (2006), C, which does not feel any radiation pressurearound β Pic, is acting as a braking agent in the form of C + via Coulomb interactions. In order to test whether theamount of carbon needed to brake metals is consistent with the new ALMA C I data, we consider an accretion diskmodel extending from 10 to 50 au with a carbon surface density Σ C ∝ r − and temperature T ∝ r − . (with T = 75 Kat 85 au) (Lynden-Bell & Pringle 1974). We compute the emission from this model as described in section 4.2. Wethen search for C I emission only in those regions of the datacube where the model predicts emission. However, wewant to exclude regions of the datacube that can contain emission from the outer disk. Thus, for those data pointswith | v | < v orb (50au) (with v orb ( r ) the orbital speed at radius r ), we request points to be at least one spatial resolutionelement inside of the line in the pv diagram defining a thin ring at 50 au (see figure 2). For points with a larger | v | , wecan be sure that no contamination from the outer disk is present. We also exclude points with a predicted emissionbelow 10% of the model peak to avoid considering regions of the datacube where only weak emission is expected anyway.We can then measure the flux in this region of the datacube, with the error estimated with a method analogous to whatis described in appendix A. Defining a detection threshold of 3 σ , we do not detect significant emission. We derive anupper limit on the C density by adjusting the model such that the probability for a non-detection, i.e. measuring a fluxbelow the detection threshold, would only be 1% if the model was correct. Assuming Gaussian noise, the probabilityof a model with flux F mod to remain undetected is given by (cid:104) (cid:16) nσ − F mod √ σ (cid:17)(cid:105) with σ the error on the flux and n = 3 in our case. The upper limit model has a mid-plane C + number density of ∼
600 cm − at 30 au, while ∼
100 cm − are necessary to brake the metals (Fern´andez et al. 2006). Thus, the data are consistent with enough carbon beingpresent in the inner disk to brake metals. However, to get better constrains, we also look at the C II observations from9 Herschel /HIFI by Cataldi et al. (2014). We measure the flux in the wings of the C II spectrum corresponding to radialvelocities between v orb (50au) and v orb (10au). Errorbars are calculated by calculating the flux in spectral regions ofthe same size without line emission and taking the standard deviation. We first detemine the errors on the flux in theleft wing and right wing individually and then add them in quadrature to obtain the error on the total measured flux.No significant (larger than 3 σ ) flux is detected neither in the H nor V beam of the data. Thus, we adjust the modelsuch that the combined probability to remain undetected in the ALMA C I and the HIFI C II (H and V beam) data isonly 1% (including the ALMA data has a negligible effect as the HIFI data are more constraining). This model has amid-plane C + density at 30 au of ∼
380 cm − . We conclude that the currently available data are consistent with carbonbeing present inside of 50 au at a level that is sufficient to brake metals. We suggest that this gas in the inner regionwas not produced in the same event as the gas seen at larger distances. If it was, we would expect a full accretion disk,which is not supported by the data (see also Cataldi et al. 2014). Instead, it might be the leftover from a previousgas-producing event. 5.5. Origin of the observed C asymmetry
A key result from our new ALMA data is that C shows the same asymmetry as CO: a clump on the SW side of thedisk. This is surprising since one would expect C to spread in azimuth within a few orbits even though C productionmight happen preferably at the CO clump. This is in contrast to CO which has a lifetime shorter than an orbitalperiod and thus remains asymmetric. Here we discuss possible solutions to this puzzle.5.5.1.
Recent event
Perhaps the simplest explanation for the observed C asymmetry is that C production at the clump location startedso recently that there was not yet enough time for azimuthal spreading (respectively symmetrisation). We expect thetimescale for azimuthal symmetrisation to be on the order of a few orbits (at 85 au, the orbital period is ∼
600 years).We investigate the symmetrisation with a 1D toy model, where the only dimension is the azimuthal angle φ . To start,we assume that all C is produced in a single point at azimuth φ , at a distance of 85 au, with a rate equal to theCO destruction rate calculated in section 5.3. In reality, only 30% of the CO emission is found in the clump (Dentet al. 2014). This setup thus maximises the asymmetry, so the derived symmetrisation timescale can be considered anupper limit. We then write the following simple equations describing the temporal evolution of the neutral and ionisedcarbon densities under the influence of C production, ionisation, recombination and orbital motion: ∂ρ ∂t = δ ( φ − φ )Λ CO m C m CO + ρ + n e γ − ρ Γ − ω ∂ρ ∂φ (5) ∂ρ + ∂t = − ρ + n e γ + ρ Γ − ω ∂ρ + ∂φ (6)where ρ and ρ + are the azimuthal densities (in kg rad − ) of neutral and ionised C respectively, Λ CO is the COdestruction rate (in kg s − ), m C and m CO are atomic and molecular masses respectively, n e the electron density, γ the recombination coefficient, Γ the ionisation rate and ω the angular velocity at 85 au. For the electron density, weassume that all electrons come from the ionisation of C, and that the typical volume extends over ∆ r =70 au in theradial direction (best-fit ‘ring’ model, section 4.2.1) and over ∆ z = 2 H (with H the scale height, see equation 2) in thevertical direction (the model remains one-dimensional; we only assume a certain volume to get a value for the electrondensity). This means that n e = ρ + / ( m C + r ∆ r ∆ z ) where m C + is the mass of a C + ion. Then, equations 5 and 6 areused to compute the change of the azimuthal densities at each time step. Figure 13 shows how the SW/NE mass ratioof neutral carbon evolves with time. The wavy pattern is due to the orbital motion of the gas. The local minimacorrespond to the times when the first gas produced by the point source leaves the NE side and enters the SW side.After ∼ years, the ratio falls below the value of our best-fit ‘ring + clump’ model.We can also use this model to estimate the time it will take to symmetrise the C from the state that what weobserve today if there is no mechanism that prevents symmetrisation. For example, assuming that only 30% of the Cproduction occurs at the clump, with the other 70% uniformly distributed at all azimuths (this corresponds to addinga production term independent of φ to equation 5), we get the green curve shown in figure 13. As expected, thesymmetrisation occurs faster—within ∼ years, the SW/NE mass ratio falls below 1.1.A caveat to this toy model is the possibility that the gas production is not constant over time. Also, for more robustconstraints, detailed hydrodynamical calculations would be necessary.0 t [yr]1234567 n e u t r a l c a r b o n S W / N E m a ss r a t i o point source30% point source, 70% uniform sourcering + clump model Figure 13.
Temporal evolution of the NE/SW ratio of the neutral carbon mass in our toy model. Two scenarios are considered:production from a point source (blue line), providing an upper limit on the symmetrisation timescale, and a scenario with initialconditions similar to what we observe today (green line). The dashed line denotes the mass ratio of our best-fit ‘ring+clump’model.
The upper limit on the symmetrisation timescale is below the lower bound of the timescale range over which Cproduction occurred, as derived from the C mass (section 5.3). Thus, it appears unlikely that the asymmetry can beexplained by invoking a very recent event.In summary, the estimated lifetime ( ∼ Resonance trapping by planet
Matr`a et al. (2017b) argued that the asymmetry in β Pic, and perhaps in a number of other debris disks, is theresult of planetesimals being trapped into mean-motion resonances (MMR) by a (migrating) planet. For β Pic, ourresolved C I observation excludes this possibility.For planetesimals in resonance with a planet there would be an enhancement in particle density at some specialazimuth relative to the planet. In the case of 2:1 MMR, such a resonance island lies 90 ◦ behind and ahead the azimuthof the perturbing planet. As collisions are more frequent in the denser region, it is expected, in this scenario, thatCO, which essentially traces recent collisions because of its short lifetime, is enhanced near the island and thereforeasymmetric. However, as the planet revolves in its orbit, the resonance island sweeps through all azimuths. Integratedover many periods, the resonance island does not linger particularly long around a special azimuth. So if we look ata tracer gas that is the integrated production of CO over many periods, as the carbon gas is (produced over ∼ ∼
10 orbits), we expect it to be evenly spread out in the orbit. The observed asymmetry in C thus cannot beexplained by its parent planetesimals being trapped into a planet MMR.Can carbon gas also be trapped into a MMR? In contrast to dust grains, gas cannot be trapped into an MMR. Onecan show that the resonance width, even when forced by a highly eccentric Jovian planet, is narrower than the vertical1scale height of the gas ( ∼ a few au). As the vertical scale height corresponds to the sound speed travel time over anorbit, gas pressure can, in an orbital period, easily disperse the gas azimuthally and radially over an extent wider thanthe resonance width. It is difficult for the weak planetary perturbation to restrain them.In section 4.2.3 we showed that an eccentric gas distribution is qualitatively able to reproduce the carbon observations.Thus, in the following sections, we discuss how such an eccentric disk could arise.5.5.3. Initially eccentric disk
If the parent planetesimal disk is eccentric because it has been left in that state by some unknown initial condition,the debris disk will initially be asymmetric. Such a disk, if precessed rigidly, can retain its initial configuration overa time much longer than the sound crossing time. However, over the lifetime of the system (20 Ma), differentialprecession should have disturbed this initial imprint. Ignoring the presence of other planets, and only considering thequadrupole precessional effect of β Pic b (Lagrange et al. 2009), the timescale for order unity change in the relativeprecession angle is (Murray & Dermott 2000) t smear ∼ M ∗ M p (cid:18) aa p (cid:19) (cid:18) aa (cid:19) − P orb ∼
14 Ma , (7)where we evaluate the expression using M p = 13 M Jup (Morzinski et al. 2015) for the planet mass and a p = 9 au(Wang et al. 2016) for its semi-major axis with the ring at a ∼
85 au and with a fiducial radial width of ∆ a ∼
20 au(consistent with that in Table 2). So even without an additional planet, an initially eccentric disk is expected to belargely smeared out. 5.5.4.
Eccentric disk secularly forced by a planet
Secular perturbation from a planet can not be responsible for producing the eccentric gas disk that we observe. Theage of the C gas ( ∼ ≥ P orb /µ , where µ is the mass ratio ofplanet to star). So any eccentricity in the gas disk will have to be inherited from their planetesimal parent bodies.But if so, what could possibly have started the parent bodies’ grind-down a mere 5 000 a ago, if they have lived in suchstates for a secular timescale? But let us ignore this issue here and proceed to consider the geometry.In the hypothetical case of long-lived parent bodies forced to a relatively low eccentricity, more particles will be foundnear apoaps where they move the slowest (the so-called ‘apo-centre’ glow, Pan et al. 2016), while collisions are morelikely to occur near the pericentre where the particle streamlines are more densely spaced and their relative velocitiesare higher. Collisions near the pericentre occur at a higher relative velocity, allowing smaller debris (which are morepopulous and have a larger collisional surface area) to break apart a given fragment. Matr`a et al. (2017a) suggestthat the mass loss rate is enhanced at either periaps or apoaps depending on the proper eccentricity and strengthof the planetesimals. From preliminary numerical simulations, we favour an enhancement at periaps (more detailedsimulations will be presented in a forthcoming paper). As a result, we expect CO (which reflects the instantaneouscollisional mass loss rate) to be concentrated in periaps (the SW side, as is observed), while both the submm emissionand the carbon line fluxes should be determined mostly by particle trajectories and should peak at apoaps, contraryto the observations (the observed C I periapsis to apoapsis flux ratio is 1 . ± .
2, see section 3.1). This argument iseasy to understand if one thinks of the CI gas as exact analogue of the small dust grains, which are also integratedproducts of previous collisions. Small grains will shine more brightly in apo-centre because there are more of themthere—unless their eccentricities are so large that the fall-off of stellar-flux at apo-centre is more important than thenumber excess, see our discussion below.At higher eccentricities, the apocentre is much further than the pericentre, and the drop-off in stellar flux, particlevolume density and gas temperature could compensate for the above trajectory effect and reduce the apocentre glow.Colder submm particles have reduced emissivity. The lower gas temperature and volume density also mean that thecarbon gas is less excited, leading to reduced line fluxes. This effect is more severe when the orbits have an intrinsicspread in eccentricity, leading to a spread in apocentre distances and further reducing the dust and gas volume densitiesin those locations. Carbon ionization fraction, on the other hand, may also be affected by the lower ionizing flux andthe smaller recombination rate.So to explain the observed same-sided asymmetry in different tracers, we will require a medium to high eccentricity( e (cid:38) .
3, table 2) and preferably a large spread in eccentricities (discussed below). For a secularly forced ring, theforced eccentricity e ∼ / a p /a ring ) e p , where e p is the eccentricity of the planet (Murray & Dermott 2000). A planetwith considerable eccentricity is then required, which may itself decimate the planetesimal belt by dynamical ejection.2Another argument against the planet hypothesis comes from the range of eccentricities required to explain theobservations. For a disk with a 20% spread in semi-major axis (∆ a ∼
20 au if a = 85 au), we expect a 20% spreadin the forced eccentricities. This is too small to help explain the same-side asymmetry and it is also smaller thanindicated by our best-fit solution.In conclusion, it seems the hypothesis that an underlying highly-eccentric planetesimal belt, forced secularly by aplanet, is unlikely to explain our observations. 5.5.5. Giant Impact
Jackson et al. (2014) proposed a scenario where a giant impact between two comparably massed bodies producesa wide spread of debris that have a range of eccentricities but a similar alignment. This is qualitatively plausible toexplain the observations (but see Matr`a et al. 2017b, for a counter argument regarding the radial width). If such anevent occurred in the recent past, it provides an interesting explanation for our deduced carbon production timescale.However, such an event seems exceedingly improbable, as we will show with an order-of-magnitude estimates of theevent rate.Consider N bodies with radius R , in a belt of semi-major axis a and width ∆ a , performing vertical epi-cyclesaround the mid-plane. Viewed by each one of these bodies, the other bodies present a certain optical depth of τ ⊥ = ( N πR ) / (2 πa ∆ a ), or a mean collision time of P orb /τ ⊥ . Summed over N bodies, this implies a mean event timeof t coll ∼ P orb πa ∆ aN × πR ∼ × yrs × (cid:16) a
80 au (cid:17) (cid:18) ∆ a
20 au (cid:19) (cid:18) N (cid:19) − (cid:18) R (cid:19) − , (8)where we assume there is no gravitational focusing among these low-mass objects, reasonable if their velocity dispersionlies much beyond their surface escape velocities. The scaled value R = 2000km is appropriate to explain the totalgas/dust mass observed, and ∆ a = 20 AU is inspired by the best fit in Table 2. The value N = 1000 is a place-holder.So, to produce an event as recent as 5 000 a ago, one would require N ∼ , or an absurdly high total mass of ∼ M ⊕ . This argument makes giant impacts exceedingly implausible at this location in the disk. Another difficultywith such a scenario is that giant impacts tend to be completely accretionary at low relative speeds, and even atspeeds beyond the surface escape, only a small fraction of mass can be unbound and released into the circumstellarenvironment (Agnor & Asphaug 2004). 5.5.6. Tidal Disruption
Here, we briefly propose an alternative scenario to produce the observed disk. A more detailed calculation will bepresented in an upcoming paper. Consider that the outer disk contains a number of Neptune-like planets ( N N ∼ afew). They have a surface escape velocity of v esc ∼
23 km s − . Let us also assume that there are N bodies similarin mass to the Moon ( ∼ . M ⊕ ) or Mars ( ∼ . M ⊕ ) and moving in space with a relatively low dispersion velocity, σ (cid:28) v esc . There is little gravitational focussing among these bodies, but strong focussing by the Neptunes. Thetimescale for a physical collision with a Neptune is t coll ∼ P orb πa ∆ aN N N × πR N (cid:18) σv esc (cid:19) ∼ × yrs × (cid:16) a
80 au (cid:17) (cid:18) ∆ a
20 au (cid:19) (cid:18) N (cid:19) − (cid:18) σ − (cid:19) (cid:18) v esc
23 km s − (cid:19) − . (9)where we chose N N = 5. This is somewhat arbitrary, but the abundance of cold Neptunes is already suggested bymicro-lensing studies, which claim that their abundance rate per star (mostly M-dwarfs) is 52% (Cassan et al. 2012).This timescale is becoming astrophysically interesting. But a physical collision with a Neptune will not produce adebris disk around the star. Instead, we focus on encounters that are close enough that the body is tidally disrupted.This means encounter distance of ∼ R N and the encounter time t disruption ∼ t coll / ∼
10 Ma (at 2 R N , while thegeometric cross-section goes up by a factor of 4, the gravitational focussing factor, ( σ/v esc ) , evaluated at R = 2 R N goes down by a factor of 2). So there could have been a few tidal disruption events over the lifetime of β Pic, if thereare thousands of moons floating around and if these moons retain low enough dispersion velocities to experience stronggravitational focussing by Neptune. This is still well above the 5 000 a event time we infer and still presents a tension.The end product of the tidal disruption shares many similarities with that from a giant impact. First, the debris willbe ejected on a variety of orbits, bringing about a large spread in eccentricity. The semi-major axis will also have aspread. All debris will return to the disruption site, making this location a region of frequent collisions. The size of thedisruption site, however, is larger than that in giant impact. As the debris flies away from Neptune, its stellar-centric3orbit is altered by Neptune’s gravity while it is still within Neptune’s Hill sphere. As a result, unlike the narrow nozzleof the size of the impactor radius in the case of giant impacts, here, the nozzle has a typical spread of order the planet’sHill radius, which reduces the peak collision rate. However, Matr`a et al. (2017b) measured a radial extent of the COclump of at least 100 au, which is clearly larger than the planet’s Hill sphere. More detailed modelling is needed tosee whether a tidal disruption event can produce such a radially broad clump. Also, the clump in the deprojection ofMatr`a et al. (2017b) might appear more extended than it is in reality because the intrinsic velocity dispersion of thegas and the finite velocity resolution of the instrument degrade the resolution of the deprojection in the y direction.In addition, the deprojection is carried out assuming velocities corresponding to circular orbits. Thus, if the orbits areeccentric in reality, the deprojection will be distorted and not correspond to the actual gas distribution.The event in β Pic is recent and must also be relatively short-lived, assuming we are not observing it at a specialtime. The duration of a tidal disruption flare will depend on the above-discussed collisional geometry, but also on thedistribution of the largest fragments (that remain or reform) after the disruption event. Both of these effects need tobe studied in detail. In addition, a short lifetime will also ensure that debris products from previous events do notinterfere with the current event. If not, the asymmetry would be washed out by previous debris which likely have adifferent asymmetry; and we would observe the radially diffused accretion disk from gas produced in previous events.In summary, our observation disfavours a few proposed scenarios (planet MMR trapping, giant impact, secularforcing). Instead, we propose that the bright, asymmetric debris disk in β Pic could be the result of a recent tidaldisruption of a Moon to Mars-sized object by a Neptune-like planet. Using the C mass derived in section 4.1, andassuming that the disrupted body has had a CO mass fraction of 10%, its total mass would indeed be (cid:38) M Moon (lowerlimit because not all carbon might have been released yet).5.6.
Comparison with the CO emission
Since C and CO have similar horizontal emission profiles (see Figure 3), the question arises how similar the spatialand spectral distribution of the emission really is. To answer this question, we first interpolate the CO data cubeonto the coordinates of the C I data cube. Then, we use convolution to adjust the spatio-spectral resolution of thedata cubes. Finally, we normalise the data cubes and subtract. Figure 14 shows the moment 0 and pv diagram ofthe residual cube for the CO 2–1 and 3–2 transitions (Matr`a et al. 2017b). The CO emission seems similar to the C I emission, with few residuals above 3 σ seen. With the data at hand, we are thus not able to detect clear differences inthe emission distribution. Note that there is significant difference in the distribution of CO 2–1 vs CO 3–2 emission(Matr`a et al. 2017b). 5.7. Comparison to the spatial distribution of metals
Brandeker et al. (2004) and Nilsson et al. (2012) found that Na and Fe have a NE/SW asymmetry reminiscent ofwhat we see for C and CO (compare figure 3 of this work to figure 3 in Brandeker et al. 2004). Na and Fe also showa tilt similar to what is seen for C and CO. The radial distribution is quite different, however, with the density beinghigher closer to the star instead of peaking at 85 au. The asymmetry is seen in both the inner and outer parts of thedisk; concerning the inner distribution of Na and Fe, the emission is traced inwards to the observational limit of 13 auin the NE, while to the SW the density peak appears to be located much further out, beyond 50 au. In the outer partsof the disk, Na and Fe can be seen all the way to the limit of the observation at the projected distance of 330 au inthe NE, while the disk seems truncated in the SW at 150–200 au.A possible scenario is that the origin of CO (and thus C and O as its dissociation products) is different from theorigin of the metals observed in the optical (Na, Fe, Ca, Ni, Ti, and Cr, Brandeker et al. 2004). While CO likely comesfrom the disruption of CO-rich cometary bodies at 85 au, the metals observed in the optical could be produced by theso-called falling evaporating bodies (FEBs, e.g. Vidal-Madjar et al. 2017) close to the star and then diffused outwardsby the stellar radiation pressure. As the presence of C + would act as a braking agent (Fern´andez et al. 2006), its spatialdistribution would naturally become imprinted on the spatial distribution of the metals. The spatial distribution ofNa and Fe is therefore consistent with C + being in an eccentric distribution resulting from a tidal disruption event,as outlined in sections 4.2.3 and 5.5.6. With the SW clump at 85 au being at the convergence point for a family ofeccentric orbits, the C would be more concentrated to the clump location in the SW, while being much more radiallydistributed in the NE, in agreement with the observed Na and Fe distributions. SUMMARY AND CONCLUSIONS4
Figure 14.
Residual moment 0 and pv diagram when subtracting the peak-scaled CO from the C I data cube. The spectro-spatial resolution has been adjust prior to subtraction. Contours are drawn at intervals of 3 σ . The residual data cube wasintegrated over ± − for the moment 0 and ±
30 au for the pv diagram respectively.
In this paper, we present resolved ALMA band 8 observations of C I emission towards the β Pic debris disk. Ourwork can be summarised as follows:1. Using a simple 1D model that calculates the ionisation balance and non-LTE level populations, we estimate thetotal C mass to be between 5 × − and 3 . × − M ⊕ . Assuming that C is produced from the photodissociationof CO at a constant rate, and that C is not removed from the system, this mass implies that gas productionstarted only ∼ I shows the same asymmetry as seen for CO: a clump on the SW side of the disk. By modellingthe spatial distribution of the C gas, we find that a satisfactory fit to the C I and archival C II data can befound by assuming that the gas consists of a ring between 50 and 120 au with a superimposed clump at the samelocation as the CO clump. A model assuming eccentric orbits of the gas with a flat eccentricity distributionbetween 0 and 0.3 also reasonably fits the data.3. The C I data are not consistent with the accretion disk predicted by Kral et al. (2016). If C gas productionindeed started only ∼ ∼ β Pic is eccentric and might have originated from a recent tidal disruption event. This could potentially alsoexplain the asymmetry observed in Na I and Fe I by Brandeker et al. (2004).A detailed study of the tidal disruption mechanism will be presented in a forthcoming paper. Follow-up observationsof the C I line at higher SNR and spectro-spatial resolution can be used to confirm or reject our hypothesis of aneccentric disk due to a tidal disruption event.We are grateful to the anonymous referee for his/her exceptionally thorough review, leading to a significant improve-ment of our manuscript. We would like to thank S´ebastien Muller and Ivan Mart´ı-Vidal from the Nordic ARC nodefor their crucial support in calibrating the data. We acknowledge helpful discussions with Thayne Currie, QuentinKral, Luca Matr`a, Nagayoshi Ohashi, Alfred Vidal-Madjar and Mark Wyatt. Luca Matr`a kindly provided the ALMACO 3–2 and 2–1 data as well as the model of the dust radiation field. Quentin Kral helped getting the LIME code run-ning. This work was partially supported by JSPS KAKENHI Grant Number JP16F16770. This work was supportedby the Momentum grant of the MTA CSFK Lend¨ulet Disk Research Group. This paper makes use of the followingALMA data: ADS/JAO.ALMA
Facilities:
ALMA, Herschel (PACS, HIFI)
Software: astropy (Astropy Collaboration et al. 2013), CASA (McMullin et al. 2007), emcee (Foreman-Mackey et al.2013), LIME (Brinch & Hogerheijde 2010), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), pythonradex( https://github.com/gica3618/pythonradex ), SciPy ( )APPENDIX A. DETAILS OF THE ERROR CALCULATIONIn this section we describe how the errors quoted in section 3.1 and 3.3 are derived.The total line emission was measured by integrating the data cube within ± − in the spectral dimesion andover the rectangular box shown in Figure 1 in the spatial dimension. Having measured the noise in the data cube, anaive way to estimate the error σ F on the total flux would be to write σ F = √ nσ dΩd v with n the number of datapoints over which the integration extends, σ the noise in the data cube and dΩ and d v the extent in solid angle andvelocity of a single data point. However, since neighboring data points are correlated, this approach is not valid.Instead, we consider the flux of a number of volumes with the same size as the volume used to measure the total flux,placed in regions of the data cube without emission. Taking the standard deviation of these flux samples, we obtain anestimate of the error. The volumes are placed such that they are sufficiently distant from each other to be consideredindependent. Since no primary beam correction has been applied, the noise can be assumed uniform over the datacube.For the continuum, the above procedure yields too few flux samples. Thus, we reduce the size of the flux samples inthe x direction to get more samples. Then, we set σ F = σ s √ N where σ s is the standard deviation of the flux samplesand N is the number of flux samples that fit into the region for which the flux is measured. Again, the flux samplesare placed sufficiently distant from each other to be considered independent.For the emission profile along the disk (Figures 1 and 3 right), we consider sample ‘volumes’ that extend only onepixel in the x direction. In the case of the profile derived from the restricted region in pv space (Figure 3), the size ofthe region over which we integrate depends on x . Thus, we take flux samples for each x individually. For this profile,the error thus depends on x .6 B. MEASUREMENT OF THE LINE BROADENING PARAMETERThe line broadening parameter b (defined as b = FWHM / (2 √ ln 2) where FWHM is the full width at half maximumof the line) parametrises the line broadening due to effects other than the orbital motion of the gas (e.g. thermalbroadening or turbulence). We can use the ALMA observations of C and CO to measure b by considering a verticalcut in the PV diagram (Figure 2) going through x = 0, i.e. considering the line of sight towards the star. For circularorbits, all gas along this line of sight is centred at the same radial velocity (namely 0 km s − ), allowing a measurementof b . However, if orbits are eccentric, the line of sight towards the star can contain additional broadening due to orbitalmotion. Furthermore, the finite spatial resolution of the instrument will blur material with non-zero radial velocityinto the line of sight towards the star. Thus, our derived value of b should be considered an upper limit.We consider the PV diagrams of C I (this work), CO 2–1 and CO 3–2 (Matr`a et al. 2017b) and for each of them fita Gaussian to the vertical cut through x = 0, yielding three independent measurements of b = (cid:112) b − b where b fit is the broadening parameter of the fitted Gaussian and b inst is the broadening due to the spectral response functionof the instrument . We estimate errors using emcee , a python implementation of an affine invariant Markov chainMonte Carlo (MCMC) sampler (Foreman-Mackey et al. 2013). The likelihood was defined via a χ measure analogousto equations 1 and 2 by Booth et al. (2017). In particular, we also use a noise correlation ratio, defined in our case asthe square root of the number of spectral pixels per spectral resolution element. We consider invariant, uninformativepriors, imposing only that the peak of the Gaussian and the broadening parameter are positive. We use 100 walkerswith 10 steps. This yields the following results for b : 0 . ± .
11 km s − (C I ), 0 . ± .
04 km s − (CO 2–1) and0 . ± .
05 km s − (CO 3–2). Combining these measurements as a weighted mean yields b = 0 . ± .
03 km s − .This value may be compared to previous measurements of b . Jolly et al. (1998) measured 0 . ± .
05 km s − from COabsorption (they also measured a high b value of 4.2 km s − from C I absorption, potentially caused by modellingdifficulties because the line is saturated). Roberge et al. (2000) obtained 1 . ± . − from C I absorption and1 . ± . − from CO absorption. Cataldi et al. (2014) and Nilsson et al. (2012) previously adopted 1.5 km s − for their models based on measurements of Ca II K absorption (Crawford et al. 1994). The b value derived from theALMA data seems thus generally lower than in previous publications. C. DETAILS OF THE DEPROJECTION PROCEDUREThe pv diagram can be used to get a deprojected view of the emission in the ( x, y ) plane. As is discussed in appendixC of Matr`a et al. (2017b), a radius r = ( GM ∗ x /v ) / can be assigned to points ( x, v ) of the pv diagram for an edge-ondisk and circular Keplerian orbits. From this, we can find the position along the line of sight y = ±√ r − x . Notethat for some points ( x, v ) of the pv diagram, the term under the square root becomes negative. This simply meansthat for the given x , the radial velocity v cannot be reached anywhere along the line of sight, i.e. we have | v | > v max ,where v max is the orbital velocity of the orbit with r = x . One has to choose how to distribute the flux from a givenpv point among the two possible y points (in front and behind the sky plane).In practice, it is easiest to take an inverse approach. First, we decide to place all flux in front of the sky plane (i.e.we only consider y < x, y ), we calculate r = (cid:112) x + y . From this, v = − (cid:113) GM ∗ r xr , wherethe minus sign accounts for the known rotation sense of the β Pic disk. We then look up the flux at ( x, v ) in the pvdiagram and assign it to the point ( x, y ) in the deprojection.In order to get a deprojection with a straightforward interpretation in terms of the distribution of the flux, it isalso necessary to transform the flux units from W m − Hz − rad − (as in the pv diagram, see figure 2) to W m − sr − .To do this, it is helpful to see the deprojection as a coordinate transformation from ( x, v ) to ( x, y ). The Jacobiandeterminant of this transformation (for y = −√ r − x as in figure 4) reads J = 32 (cid:112) GM ∗ xy ( x + y ) − / (C1)The flux read from the pv diagram is then multiplied by | J | and an additional constant factor ν c d (with ν the centralfrequency of the line and d the distance between the observer and β Pic).If the flux units of the deprojection are not transformed (as for example in Matr`a et al. 2017b), an interpretationof the deprojection is not straightforward because the flux in a certain area of the deprojection does not equal the See ALMA Cycle 6 Technical Handbook, section 5.5.2 (Spectral Setup), https://almascience.nao.ac.jp/documents-and-tools/cycle6/alma-technical-handbook Figure 15.
Same as figure 4, but in the same units as the pv diagram, i.e. without multiplying the deprojection by the Jacobianshown in equation C1. integral over x and y over this area. Figure 15 shows the deprojection without multiplication by the Jacobian (i.e. in thesame units as the pv diagram). Comparing to figure 4, it becomes apparent that a deprojection without transformedunits might lead to visual mis-interpretation, for example about the relative amount of flux at large radii (say beyond150 au). Indeed, from equation C1, we see that the Jacobian gets smaller at large radii. D. APPROXIMATE CALCULATION OF THE STATISTICAL EQUILIBRIUMWe solve the SE under the simplifying assumption that the photon escape fraction is high, i.e. we neglect(de)excitation by emitted line photons and assume that the background radiation field (CMB, stellar field and dustcontinuum) is not attenuated by the gas. In this section, we justify these assumptions.The SE equation for an atomic level i can be written (cid:88) j (cid:54) = i (cid:0) x j ( A ji + C ji + ¯ JB ji ) − x i ( A ij + C ij + ¯ JB ij ) (cid:1) = 0 (D2)where x j is the fractional population of level j , C ij = K ij n e is the collisional (de)excitation rate (with K ij thecollisional rate coefficient and n e the electron density), ¯ J the frequency-integrated mean intensity and A ij and B ij arethe Einstein coefficients (where we set A ij = 0 if i < j ). If we can show that A ji + C ji (cid:29) ¯ JB ji (for any i, j ), then wehave demonstrated that background radiation and line photons are not important to solve the SE.The frequency-integrated mean intensity is the sum of the line emission and background radiation at each point ofthe disk: ¯ J = ¯ J line + ¯ J backg . To get insight into the individual contribution of the two components, we consider themseparately. First, we calculate R backg = ¯ J backg B ji / ( A ji + C ji ) for all transitions and all locations (except where thegas density is below 5% of the peak density) for the best-fitting models described in sections 4.2.1, 4.2.2 and 4.2.3. Wefind that R backg < .
02 for C I and R backg < .
001 for C II . Thus, although we included background radiation in our8calculation, it is actually negligible. We double-check by re-computing the models without background radiation andfind that the results indeed do not change.Next, we consider (de)excitation by line emission and compute R line . To this end, we compute, for every location(again, except where the gas density is below 5% of the peak density), the number of line photons arriving from theother grid points, neglecting optical depth and the velocity field (the velocity field could red/blue-shift photons fromother grid points out of the transition). We are thus calculating an upper limit on ¯ J line . We find that R line < . I and R line < .
03 for C II . Thus, a more sophisticated model should include the full radiative transfer for C I , butneglecting the line photons is a decent approximation given the quality of our data, since it considerably simplifies thecalculation of the models.As an additional test, we used the LIME code version 1.9.1 (Brinch & Hogerheijde 2010) to calculate the full non-LTEradiative transfer (neglecting background radiation except for the CMB) for our best-fit models from sections 4.2.1,4.2.2 and 4.2.3. We find that the total flux computed by
LIME differs at most by a factor 1.2 for both C I and C II . E. GAS DENSITY FOR ECCENTRIC ORBITSIn this section, we derive the gas surface density of the model with eccentric orbits presented in section 4.2.3. Weassume that all orbits have a common pericentre and that the distribution of eccentricities is known and given by P ( e ). We search the probability P ( r, θ ) (which we assume is proportional to the surface density) to find a particle atradius r and true anomaly θ . We first consider P ( e, θ ) = P ( e ) P ( θ | e ). In our model, a given eccentricity e correspondsto a single orbit (because all orbits have a common pericentre). Thus, P ( θ | e ) is interpreted as the probability to finda particle at true anomaly θ along the orbit with eccentricity e . The time a particle spends at a given θ is inverselyproportional to the orbital velocity. Thus, we have P ( θ | e ) = C ( e ) | v ( e, θ ) | (E3)where the normalisation constant C ( e ) = (cid:16)(cid:82) | v ( e,θ (cid:48) ) | d θ (cid:48) (cid:17) − . The orbital velocity is given by v ( e, θ ) = (cid:115) µ e cos θ + e q (1 + e ) (E4)with µ = GM ∗ and q the fixed pericentre distance. Next, we consider the transformation from ( e, θ ) to ( r, θ ) where r = q (1 + e )1 + e cos θ (E5)The transformation is bijective, except for θ = 0. Therefore, we can write P ( r, θ ) = P ( e, θ ) | J | (E6)with J the Jacobian determinant of the transformation, given by J = ∂e∂r , which can be calculated from equation E5.At θ = 0, the density P ( r, θ ) diverges. In practice, this is easily handled by simply not sampling the singularity on thenumerical grid. REFERENCES Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12Agnor, C., & Asphaug, E. 2004, ApJL, 613, L157Apai, D., Schneider, G., Grady, C. A., et al. 2015, ApJ,800, 136Armitage, P. J. 2009, Astrophysics of Planet Formation(Cambridge University Press),doi:10.1017/CBO9780511802225 Artymowicz, P. 1997, Annual Review of Earth andPlanetary Sciences, 25, 175Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33Aumann, H. H. 1985, PASP, 97, 885Backman, D. E., & Paresce, F. 1993, in Protostars andPlanets III, ed. E. H. Levy & J. I. Lunine, 1253–1304 Booth, M., Dent, W. R. F., Jord´an, A., et al. 2017,MNRAS, 469, 3200Bouret, J.-C., Deleuil, M., Lanz, T., et al. 2002, A&A, 390,1049Brandeker, A. 2011, ApJ, 729, 122Brandeker, A., Liseau, R., Olofsson, G., & Fridlund, M.2004, A&A, 413, 681Brandeker, A., Cataldi, G., Olofsson, G., et al. 2016, A&A,591, A27Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature,481, 167Cataldi, G., Brandeker, A., Olofsson, G., et al. 2014, A&A,563, A66Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153,182Crawford, I. A., Spyromilio, J., Barlow, M. J., Diego, F., &Lagrange, A. M. 1994, MNRAS, 266, L65Crifo, F., Vidal-Madjar, A., Lallement, R., Ferlet, R., &Gerbaldi, M. 1997, A&A, 320, L29Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014,Science, 343, 1490Draine, B. T. 1978, ApJS, 36, 595Fern´andez, R., Brandeker, A., & Wu, Y. 2006, ApJ, 643,509Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,J. 2013, PASP, 125, 306Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006,AJ, 132, 161Grigorieva, A., Th´ebault, P., Artymowicz, P., & Brandeker,A. 2007, A&A, 475, 755Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017,A&A, 602, A105Higuchi, A. E., Sato, A., Tsukagoshi, T., et al. 2017, ApJL,839, L14Hunter, J. D. 2007, Computing In Science & Engineering,9, 90Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D.2014, MNRAS, 440, 3757Jolly, A., McPhate, J. B., Lecavelier, A., et al. 1998, A&A,329, 1028Kiefer, F., Lecavelier des Etangs, A., Boissier, J., et al.2014, Nature, 514, 462Kral, Q., & Latter, H. 2016, MNRAS, 461, 1614Kral, Q., Matr`a, L., Wyatt, M. C., & Kennedy, G. M. 2017,MNRAS, 469, 521Kral, Q., Wyatt, M., Carswell, R. F., et al. 2016, MNRAS,461, 845Lagrange, A.-M., Gratadour, D., Chauvin, G., et al. 2009,A&A, 493, L21 Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010,Science, 329, 57Lagrange, A.-M., Boccaletti, A., Milli, J., et al. 2012, A&A,542, A40Lecavelier des Etangs, A., & Vidal-Madjar, A. 2016, A&A,588, A60Lee, L. C. 1984, ApJ, 282, 172Li, D., Telesco, C. M., & Wright, C. M. 2012, ApJ, 759, 81Lieman-Sifry, J., Hughes, A. M., Carpenter, J. M., et al.2016, ApJ, 828, 25Liseau, R. 2003, in ESA Special Publication, Vol. 539,Earths: DARWIN/TPF and the Search for ExtrasolarTerrestrial Planets, ed. M. Fridlund, T. Henning, &H. Lacoste, 135–142Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603Mamajek, E. E., & Bell, C. P. M. 2014, MNRAS, 445, 2169Matr`a, L., Wilner, D. J., ¨Oberg, K. I., et al. 2018, ApJ,853, 147Matr`a, L., MacGregor, M. A., Kalas, P., et al. 2017a, ApJ,842, 9Matr`a, L., Dent, W. R. F., Wyatt, M. C., et al. 2017b,MNRAS, 464, 1415Matthews, B. C., Krivov, A. V., Wyatt, M. C., Bryden, G.,& Eiroa, C. 2014, Protostars and Planets VI, 521McMullin, J. P., Waters, B., Schiebel, D., Young, W., &Golap, K. 2007, in Astronomical Society of the PacificConference Series, Vol. 376, Astronomical Data AnalysisSoftware and Systems XVI, ed. R. A. Shaw, F. Hill, &D. J. Bell, 127Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al.2015, ApJ, 811, 18Milli, J., Lagrange, A.-M., Mawet, D., et al. 2014, A&A,566, A91Mo´or, A., ´Abrah´am, P., Juh´asz, A., et al. 2011, ApJL, 740,L7Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015,ApJ, 815, 108Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471Murray, C. D., & Dermott, S. F. 2000, Solar SystemDynamicsNesvold, E. R., & Kuchner, M. J. 2015, ApJ, 815, 61Nilsson, R., Brandeker, A., Olofsson, G., et al. 2012, A&A,544, A134Olofsson, G., Liseau, R., & Brandeker, A. 2001, ApJL, 563,L77Pan, M., Nesvold, E. R., & Kuchner, M. J. 2016, ApJ, 832,81Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010,A&A, 518, L10