ALMA resolves C i emission from the Pictoris debris disk
The 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 CO into 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 derive a low C mass (between and M), indicating that gas production started only recently (within 5 000 years). No evidence is seen for an atomic accretion disk inwards of the CO belt, perhaps because the gas did not yet have time to spread radially. The fact that C and CO share the same asymmetry argues against a previously proposed scenario where the clump is due to an outward migrating planet trapping planetesimals in an resonance; nor can the observations be explained by an eccentric planetesimal belt secularly forced by a planet. Instead, we suggest that the dust and gas disks should be eccentric. Such a configuration, we further speculate, might be produced by a recent tidal disruption event. Assuming that the disrupted body has had a CO mass fraction of 10%, its total mass would be 3 .
0000-0002-2700-9676]Gianni Cataldi \altaffiliationOverseas researcher under Postdoctoral Fellowship of Japan Society for the Promotion of Science.
In debris disk systems, the continuous collisional destruction of larger bodies such as comets or asteroids produces abundant amounts of dust, with the smallest grains quickly removed by radiation pressure (Backman & Paresce, 1993; Wyatt, 2008). A debris disk provides evidence that planetesimal-sized bodies were formed during the earlier protoplanetary phase (Artymowicz, 1997; Matthews et al., 2014) and gives us the opportunity to study the properties of the building blocks of planets. Studying these systems is therefore intimately linked to our efforts of understanding how planets form.
The debris disk around the young ( Ma, Mamajek & Bell, 2014) A6V star (Gray et al., 2006) Pictoris has been used as a laboratory to study the early evolution of planetary systems ever since its discovery by the Infrared Astronomical Satellite (IRAS) (Aumann, 1985). Smith & Terrile (1984) obtained the first resolved image showing an edge-on disk. Since then, the properties of the dust disk have been extensively studied with observations at multiple wavelengths. Today, we know that the belt of parent bodies is located at 100 au (Dent et al., 2014) and that Pic hosts a giant planet (e.g. Chilcote et al., 2017, and references therein) with a semi-major axis of 10 au (Lecavelier des Etangs & 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 the presence of gas around Pic (Slettebak, 1975; Slettebak & Carpenter, 1983). The Pic disk is thus part of a small sub-sample of debris disks where gas has been detected. Currently, there are about 20 such gaseous debris disks known (e.g. Redfield, 2007; Moór et al., 2011; Roberge et al., 2014; Lieman-Sifry et al., 2016; Matrà et al., 2017a).
The spatial distribution of the gas around Pic has been studied with resolved observations in the optical and recently with the Atacama Large Millimeter/submillimeter Array (ALMA). These data showed that the gas disk is radially 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 absorption features 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 seen around 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 of solid material rather than leftover from the protoplanetary phase. Evidence for such a secondary scenario comes for example from theoretical arguments on the dynamical lifetime of the gas (Fernández et al., 2006). Also, models of the excitation of the CO 3–2 and 2–1 transitions observed by ALMA imply that not enough H is present in the disk to shield CO from photodissociation over the lifetime of the disk, thus the necessity of a gas replenishment mechanism (Matrà et al., 2017b). Studying this secondary gas opens up the interesting possibility to constrain the composition of the parent bodies (e.g. Kral et al., 2016; Matrà 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; Roberge et al., 2006; Brandeker et al., 2016). Recently, the first detection of hydrogen was reported by Wilson et al. (2017). CO remains the only molecule detected so-far (e.g. Dent et al., 2014; Matrà et al., 2018). The observed elemental abundances are 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ández et al. (2006) showed that the carbon overabundance provides a braking mechanism, preventing other metals that are strongly affected by radiation pressure (such as Na) from being quickly ejected from the system. Carbon also plays a crucial role in determining the excitation conditions of atomic fine-structure or molecular rotational transitions (e.g. Zagorovsky et al., 2010). This is because in a secondary, hydrogen-depleted disk, collisional excitation is likely dominated by electrons, and ionisation of carbon is the main electron source.
Using spectrally resolved observations of C ii with Herschel111Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA./HIFI, the spatial distribution of carbon was constrained by Cataldi et al. (2014). They found that most of the carbon gas is located at 100 au, with tentatively more emission from 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 with a 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 less than an orbit (Visser et al., 2009). It is thus clear that CO needs to be produced continuously, and is a source for both C and O. This might provide a natural explanation for the observed super-solar abundance of C with respect to 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. Several hypothesis have been put forward to explain the existence of the clump. Firstly, it could be the location of a giant collision. The clump then results from the fact that the orbits of the collision debris always all go through the collision point (Jackson et al., 2014). Secondly, the clump could be due to cometary bodies trapped in a resonance with an outward migrating, yet unseen giant planet (Wyatt, 2003, 2006). Using ALMA follow-up observations of the CO 2–1 transition at higher resolution, Matrà et al. (2017b) dismissed the giant collision scenario based on the large radial extent of the clump. Thirdly, Nesvold & Kuchner (2015) proposed that the interaction between a spiral density wave and 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 hydrodynamical evolution of carbon and oxygen in the disk. They assumed that the produced atomic gas is subject to the 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 first spatially resolved observations of C i, revealing its distribution in the disk. Our paper is organised as follows. We describe the observations in section 2 and present the results in section 3. In section 4, we present simple gas emission models to study the total mass and spatial distribution of the carbon gas. We discuss our results in section 5 and conclude in section 6.
2 Observations and data reduction
We observed the Pic disk using band 8 receivers of ALMA on December 19, 2015 during the ALMA cycle 3 Early Science 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) was 0.6 mm with standard deviation of 0.1 mm. The array consisted of 36 antennas arranged in a hybrid configuration containing both short and long baselines ranging from 15 m to 6.3 km. In principle, we are thus sensitive to angular scales between 0.02″ and 5″ on the sky222see ALMA Cycle 6 Technical Handbook, section 3.6 (Spatial Filtering), https://almascience.nao.ac.jp/documents-and-tools/cycle6/alma-technical-handbook. However, the visibilities from the longer baselines were affected by large atmospheric phase fluctuations and therefore flagged during the calibration process (see below). The observations were executed as a mosaic to obtain uniform sensitivity over the whole disk. One pointing was centred on the star and two additional pointings were centred at along the position angle of the disk. The size of the primary beam is 11.8″. Antenna elevations varied between 27° and 60°.
We placed three spectral windows, each with 1920 channels and a channel spacing of 488 kHz (total bandwidth 937.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 at 477.50 GHz. For C i, this corresponds to a channel spacing of 0.30 km s and an effective spectral resolution333see ALMA Cycle 6 Technical Handbook, Table 5.2 of 0.34 km s (spectral averaging factor ). In addition, a fourth spectral window with 128 channels and a total bandwidth of 1875 MHz was placed at 479 GHz in order to observe the dust continuum.
The data were calibrated using CASA 4.7.0 (McMullin et al., 2007). We performed standard water vapour radiometre (WVR) calibration and system temperature corrections. The following calibration sources were observed: J0522-3627 (bandpass), J0538-4405 (phase), J0519-4546 and J0538-4405 (flux). However, we discarded the data from the two flux calibrators and instead used J0522-3627 to calibrate both bandpass and flux because this latter source had a significantly better measurement of its absolute flux in the ALMA Calibrator Source Catalogue.
The antenna time-dependent gain calibration solutions were adversely affected by the high atmospheric phase fluctuations on the longer baselines. We therefore flagged all baselines longer than 2 km, effectively removing one third of the baselines. In addition, two bad antennas (DA44 and DA54) were also flagged. This allowed us to derive acceptable calibration solutions.
For the spectral windows placed onto emission lines, we performed continuum subtraction using the uvcontsub task within CASA. We then imaged the visibilities using the CLEAN task within CASA. In order to increase our surface brightness sensitivity, we applied a taper of 1″, thus significantly reducing the contribution of the remaining long baselines (at 492 GHz, an angular scale of 1″corresponds to a baseline length of 126 m). We also produced a continuum image using CLEAN with the same taper, combining all spectral windows except the one centred onto CS 10–9, which is 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 significantly worse than requested. Consequently, the data initially did not pass Quality Assurance 2 (QA2). However, we still decided to publish the data as the signal-to-noise ratio (SNR) is sufficient to provide new information on the Pic system.
3.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 km s (with respect to the systemic velocity). The figure has been rotated to align the horizontal direction with the main dust disk (position angle of , Lagrange et al., 2012). The beam size is 1.18″0.95″ (23 au19 au) with a major axis position angle with respect to the North of 44°. We achieve a 1 sensitivity of W m Hz sr (70 mJy beam) in the individual channels. We perform photometry by considering a rectangular aperture extending 115 au in the horizontal direction (measured from the stellar position, see Figure 1) and 30 au in the vertical direction (measured from the mid-plane). This yields a total flux of W m ( Jy km s), where the error is random without any systematic calibration uncertainty taken into account. We did not correct for the primary beam, as it changes the total flux only within the quoted error bars (the same applies for the continuum, section 3.3). The measured flux is consistent with the value of W m ( Jy km s) derived from single dish observations by Higuchi et al. (2017). The same method yields upper limits of W m (3.6 Jy km s) for CS 10–9 and 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 asymmetry with a peak on the SW side of the disk. The same asymmetry is seen for the CO emission (Dent et al., 2014; Matrà et al., 2017b) and tentatively also for C ii (Cataldi et al., 2014). However, the SW/NE flux ratio within the rectangular box of Figure 1 is not significant at (we calculate the error on the ratio by propagating the error like this: ). The SW/NE ratio of the peaks in the emission profile is . Thus, there is not necessarily more 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 direction within au). This figure thus shows the radial velocity of the emission (with respect to the systemic velocity of the star, assumed to be km s, Brandeker, 2011) as a function of the projected position along the disk mid-plane. By using the pv diagram, we can constrain the radial distribution of the emission despite the edge-on orientation of the disk. Indeed, figure 2 also shows the radial velocity for circular Keplerian orbits at 50 au and 220 au around a 1.75 M star (Crifo et al., 1997), seen edge-on. This is the approximate radial extent of the CO as found by Matrà et al. (2017b). As can be seen, no significant C i emission is detected beyond these lines, suggesting that 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, when integrating over the spectral axis, we only take data points within this region into account (rather than everything within km s as was done to produce Figure 1). The region is illustrated by the dashed lines in Figure 2 and the resulting moment 0 map and emission profile are shown in Figure 3. The SW/NE flux ratio (within the same box as in figure 1) is now . Most importantly, the SNR of the emission profile is significantly improved and the SW/NE asymmetry more clearly visible. We measure a SW/NE peak ratio of . Thus, the significance of the peak ratio asymmetry is only marginal. However, in Figure 3 (right) we also show the profiles of the CO 3–2 and 2–1 emission (see Figure 2 of Matrà et al., 2017b), for which the SW/NE flux ratios are (CO 2–1) and (CO 3–2) and the peak ratios are (CO 2–1) and (CO 3–2). The C i emission follows the CO 2–1 emission surprisingly well, with the same distinct peak on the SW side of the disk and the same asymmetric extent (out to 150 au in the NE and 100 au in the SW). This suggest that the asymmetry is indeed real. This is a surprising result. The asymmetry in CO can be readily understood from its short (less than one orbit) lifetime due to photodissociation: 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 spread in azimuth within a few orbits. In section 5.5 we discuss possible solutions to this puzzle.
From the moment 0 maps, it is also apparent that the observed emission is slightly tilted, with the NE side below and the SW side above the midplane of the main outer disk (defined by ). A similar tilt is seen for CO (Dent et al., 2014; Matrà et al., 2017b). As is discussed by Matrà et al. (2017b), two reasons might be imagined for this observation: either the gas disk is indeed tilted with respect to the dust disk, or the tilt is a projection effect, resulting from 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 inner dust 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-on view of the emission (e.g. Dent et al., 2014; Matrà et al., 2017b). Indeed, each point in the pv diagram corresponds to two points in the xy-plane of the disk (where we define as the coordinate running along the line of sight), one in front and one behind the sky plane. There remains the degeneracy of how to distribute the flux of a given pv point among the two points in the xy-plane. As discussed by Matrà et al. (2017b), the degeneracy can be broken if the disk is not perfectly edge-on. However, for simplicity, we follow Dent et al. (2014) and assume that the disk is edge-on. Our primary interest is to illustrate the radial extent of the emission and the position of the clump. We thus choose to place all flux in front of the sky plane, but other physically motivated choices are possible (Dent et al., 2014). Figure 4 shows the deprojection. Points of the pv diagram with au were masked because the radial velocity in this regions tends towards zero for all radii, i.e. it becomes difficult to assign a radius to a given radial velocity. Emission is seen approximately out to 150 au. The clump appears at a similar position angle as seen in CO (Dent et al., 2014). Details of 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 -plane but rather how much emission the observer receives from various locations in the -plane.
The continuum image at 485 GHz is shown in figure 5. The 1 noise level is W m Hz sr (1.3 mJy beam). The beam size is 1.19″0.96″ (23 au19 au) with the position angle of the major axis being 46°. We measure a total flux of W m Hz ( mJy) in the rectangular region ( au from the star along and 30 au from the mid-plane along ) indicated in the figure (see appendix A for details on the error calculation). The measured flux is consistent with a Rayleigh-Jeans extrapolation of the flux measured by ALMA at 870 m (Dent et al., 2014), and is a factor of 2 below of what is expected from the infrared to mm SED fit by Vandenbussche et al. (2010).
As was already seen in the ALMA data by Dent et al. (2014), the continuum is brighter on the SW side of the disk at projected separations between 50–100 au. However, the SW/NE integrated flux ratio is only . In any case, 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.
4.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 1D model, 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 of the 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 power beam width is only 200 au.
We need to compute the statistical equilibrium of the level populations, as in general the emission cannot be assumed to be in local thermal equilibrium (LTE). The gas in the Pic disk is thought to be of secondary origin and thus depleted in hydrogen (e.g. Matrà et al., 2017b). We thus assume that the dominant collider species is electrons, and that photoionisation of carbon is the main electron source (Fernández et al., 2006; Kral et al., 2016), i.e. the density of ionised 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 emission. The second step is to determine the amount of neutral carbon necessary to reproduce the C i flux observed by ALMA, using the electron density derived in the first step.
To solve the statistical equilibrium and radiative transfer, we use pythonradex444https://github.com/gica3618/pythonradex, a python implementation of the RADEX code (van der Tak et al., 2007) with atomic data from the LAMDA database555http://home.strw.leidenuniv.nl/~moldata/ (Schöier et al., 2005). This code uses 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 is clearly 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-major axis of 35 au in the radial direction and semi-minor axis of 10 au in the vertical direction (see Figures 1,2 and 3). Note that 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 is the dominant component. However, it turns out that including radiative (de-)excitation from these sources does not change our results because (de-)excitation is dominated by collisions and spontaneous emission.
We then consider a wide range of kinetic temperatures 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 our more detailed modelling in section 4.2, we deem lower temperatures unlikely. Detailed thermodynamical modelling by Kral et al. (2016) also suggests that the temperature is above 50 K within 100 au, although in their model, the temperature drops to 20 K at 150 au.
Figure 6 shows the derived total carbon mass as a function of . The figure also shows the individual C and C masses. To assess the importance of non-LTE effects, we also show corresponding curves for which LTE has been assumed. Both C i and C ii are generally speaking close to LTE. For higher temperatures, the C mass required to reproduce the observed C i flux is higher in LTE than in non-LTE. This simply happens because in LTE, higher levels get 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.4 for C i and 2.0 for C ii.
From Figure 6, we determine a total carbon mass between and M. The lower bound is quite robust to changes of the size of the emitting region as it is close to the LTE value. For example, increasing or decreasing the radius of the emitting sphere by 50% does not change the lower bound by more than 15%. On the other hand, it is clear that the upper bound is more uncertain as it can quickly increase if one allows for lower temperatures and/or smaller 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 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 at disposal. These more detailed models derived higher masses by using the spectrally resolved C ii observations from Herschel/HIFI: Cataldi et al. (2014) obtained M while Kral et al. (2016) derived M (the latter study also used an upper limit on the C i flux). In the Kral et al. (2016) model, the total C mass is dominated by ionised carbon that is located beyond 100 au. The C i flux is of little importance for the total C mass budget of that model. The discrepancy can thus not be explained by the fact that we include the C i measurement. On the other hand, most of the carbon in the Cataldi et al. (2014) model is located between 150 and 300 au with an ionisation fraction of roughly 50%, i.e. there is a significant contribution of neutral carbon to the total mass. However, our ALMA data show no C i emission beyond 120 au. That model indeed over-predicted the C i flux by a factor of 20, suggesting that it contains too much neutral carbon, partially explaining the difference in the mass estimates. But Cataldi et al. (2014) also derived a higher mass of ionised carbon compared to our estimate. In addition, their fit to the resolved C ii line profile (see their figure 2b) clearly suggests that C ii emission beyond 150 au is needed to fit the 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 line core). So maybe there is largely ionised carbon gas beyond 150 au present and the reason why we derive a lower ionised carbon mass here is because we assumed (based on the C i data) a too small volume (i.e. too high density and thus 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 thus not fully resolve the discrepancy, which might be due to the different assumptions of the models. The relatively small HIFI beam (FWHM of 200 au) could also play a role. For example, deviations from axisymmetry in the distribution of 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 the star.
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, OBSID 1342238190). 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 is optically thin to ionising photons. The star is the main source of ionising photons, but ionisation from the interstellar radiation field (ISRF) is also included. The ionisation balance is solved analytically by assuming that all electrons come from the photoionisation of carbon. The ionisation fraction thus only depends on the distance to the star, the local gas density and the kinetic temperature (via the recombination coefficient). Next, we use the derived electron density 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 other colliders) and radiative (de)excitation by line photons and the background radiation field, where the latter is assumed to be composed of the CMB, the star and the dust continuum. For the dust continuum, we employ the field shown in figure B1 of Kral et al. (2017) (that figure actually shows the field in the midplane of the disk at the position angle of the clump (L. Matrà, 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 approach and assume that the line emission is essentially optically thin. Basically, while the emission can become optically thick along the disk mid-plane (in our best-fit models, the maximum optical depth is typically of the order of 1), it is optically thin in the vertical direction, i.e. the photon escape fraction is high. Thus, we assume that the background field is not attenuated by the gas (all atoms are subject to the full background field) and that (de)excitation by line photons can be neglected. These assumptions make the SE easy to solve. We further discuss this approximation in appendix D. The background radiation turns out to be unimportant for our models.
Having solved the SE, we compute the emitted spectrum for each grid cell, red- or blue-shifting it according to its radial 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 models and fix 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 integrate spatially to get a model HIFI spectrum (Cataldi et al., 2014). For the C i ALMA observations, we convolve the model data cube to the same spatial and spectral resolution as the observations and multiply by the primary beam. The residual moment 0 maps (respectively pv diagrams) shown in figures 7, 8 and 11 are obtained by subtracting the model moment 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 density reads
where and are cylindrical coordinates (with perpendicular to the disk mid-plane), is the scale height, and define the radial extent of the ring and is the constant surface density. The scale height in hydrostatic equilibrium for a vertically isothermal disk is given by (e.g. Armitage, 2009)
with the Boltzmann constant, the temperature, the proton mass, the gravitational constant and the stellar mass. For the mean molecular weigh , we follow Matrà et al. (2017b) and assume . At 85 au and for 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 test different 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 correlation of neighbouring data points into account by using an appropriate noise correlation ratio for each data set). We sum the from the Herschel/HIFI data (C ii) and the ALMA data (C i).
Figure 7 shows the C i residuals in the and plane as well as the C ii line emission of the model with the lowest . The corresponding model parameters are given in table 1. The model provides a decent fit to the data, but unsurprisingly is not able to reproduce the clump observed in the SW. Thus, we consider a more complicated model by adding a clump to the uniform ring in the next section. Such a model has no direct physical justification, but is useful 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
where runs along the disk mid-plane in the sky plane and along the line of sight (and ). Furthermore, and designate the center of the clump, where we have defined as the clump’s radial distance to the star and as its azimuthal angle in the x-y-plane. We fix (Matrà et al., 2017b). The standard deviation of the clump density distribution in the --plane is and the density at the centre . The total carbon number density is then given by . We compute models over the parameter space listed in table 1. Figure 8 is the analogue of figure 7 for the best ‘ring + clump’ model. As can be seen, the model does a slightly better job in modelling the clump.
|‘ring’||‘ring + 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 trapping of cometary bodies by a migrating giant planet. However, as we discuss in section 5.5, based on our new C i data, we deem 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 of circular orbits and instead consider gas on eccentric orbits. In section 5.5.6, we discuss how such orbits could arise in the first place. As an example, we here consider a model where the eccentricity of the orbits is uniformly distributed between two values and 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 with the parameters listed in table 2. Figure 9 shows a face-on view of the mid-plane carbon gas density and figure 10 the pv diagram of the modelled C i emission. As is seen in figure 11, this model is similarly effective in reproducing the observed C i emission, although it has some difficulties to reproduce the C ii line. We emphasize that the purpose of this model is merely to demonstrate that eccentric gas orbits are able to reproduce the general morphology with the clump. 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 0.4 km s (this is smaller than and thus consistent with the broadening measured from the pv diagram, see Appendix B). The velocity shift and broadening depend on the eccentricity and orientation of the disk. To verify this velocity shift and constrain the model, we would need to know the absolute stellar radial velocity to better accuracy than 80 m s (for a 5-detection of the shift).
The column densities of neutral carbon towards the star for the three models discussed in sections 4.2.1, 4.2.2 and 4.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 cm and cm, while Roberge et al. (2006) report cm.
4.3 Absence of an atomic accretion disk
Recently, Kral et al. (2016) presented a model that computes the temperature, ionisation and hydrodynamical evolution of the atomic carbon and oxygen in the Pic disk. In this model, the atomic gas is produced in a parent belt from photodissociation of CO and then evolves viscously under the influence of the MRI, i.e. it forms an accretion and decretion disk. Kral et al. (2016) predict that the carbon is mostly neutral in the inner parts of the accretion disk.
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 the model 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 in section 4.2. Figure 12 shows the moment 0 and pv-diagrams of the model and the residuals when subtracting the model from the ALMA C i data. It is clear from this figure that the accretion profile predicted by the 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 where too much emission is present at high velocities.
5.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 (with the assumed stellar mass, Crifo et al., 1997) is preferred to reproduce the data. While this result has to 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 dynamical mass by modelling spatially and spectrally resolved Na i emission. They ascribed their finding to radiation pressure opposing the gravity of the star. In fact, radiation pressure on Na is so strong that a braking mechanism is required to keep it on the observed Keplerian orbits and preventing it from being blown out (e.g. Liseau, 2003; Brandeker et al., 2004). Fernández et al. (2006) suggested that all species affected by radiation pressure are largely ionised and coupled into a single fluid by Coulomb interactions, with carbon acting as a braking agent. In this situation, neutrals are also expected to be coupled (the only exception would be neutrals that do not ionise yet still feel a significant radiation force). Thus, one might actually expect that Na i and C i share the same dynamics. The dynamical mass (respectively the radiation pressure coefficient) is coupled to the composition of the gas (Fernández et al., 2006) and could thus give interesting information on the elemental abundances. However, given the unknown uncertainty of the derived dynamical mass which could also be influenced by the assumptions of our 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 figures 7, 8 and 11). The issue is particularly pronounced for the model with eccentric orbits. The fits to the C ii line profile by Cataldi et al. (2014) indicate that the line core requires ionised carbon beyond 150 au to be present. As already discussed in section 4.1, there might thus exist a gas component with a high ionisation fraction beyond 150 au. This component might also be required to prevent the metals observed there from being blown out by radiation pressure (see section 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 Cataldi et al. (2014) fits the C ii line profile, it would be a bad fit to our new C i data as it strongly overpredicts the C 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 carbon is 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à et al. (2017b) calculated a CO lifetime of 300 a in the clump against the ISRF. Here we revise this value, showing that photodissociation from the star actually dominates over the ISRF.
We compute the CO photodissociation rate using photodissociation cross sections from the Leiden Observatory database of ‘photodissociation and photoionisation of astrophysically relevant molecules’666http://home.strw.leidenuniv.nl/~ewine/photo/ (Visser et al., 2009; Heays et 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ández et al. (2006). This model is then complemented with data in the UV from the Space Telescope Imaging Spectrograph (STIS) aboard the Hubble Space Telescope (HST) (Roberge et al., 2000) as well as from the the Far Ultraviolet Spectroscopic Explorer (FUSE) (Bouret et al., 2002; Roberge et al., 2006). This is important because Pic shows additional emission in the UV above the predictions from a standard stellar atmosphere model. These additional UV photons impact the calculated lifetime of CO. The lifetime corresponding to the observed stellar spectrum is 70 a at 85 au. However, since the light we observe has travelled through the full CO and C column, this is actually the lifetime of a CO molecule sitting behind this 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 lifetime against stellar photons at 85 au is 20 a.
From Matrà et al. (2017b) and Roberge et al. (2000), we know the vertical column density at the clump location and the horizontal column density of CO against the star respectively. By applying a rough scaling based on the deprojected CO emission in Matrà et al. (2017b), we estimate horizontal and vertical column densities at the clump location and along the line of sight to the star. For C, the horizontal column density against the star is taken from Roberge 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 shielding towards the star, it is not difficult to see that the average CO molecule experiences half of the total column density against the star. For a Gaussian sphere, the average column density per molecule towards the ISRF turns out to be around half of the column density seen from the centre of the sphere. We choose shielding factors corresponding to these average column densities. Applying these shielding factors to the unshielded lifetimes derived above, we find that 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 M derived by Matrà et al. (2017b), a CO lifetime of 50 a leads to a CO mass loss rate of kg s. Under the assumptions that 1) the CO mass loss rate is constant, 2) C is produced only from CO photodissociation (C might also be produced from the destruction of other carbon-bearing molecules such as for example 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 in section 4.1 to calculate the time-scale over which CO and C production has been ongoing: with the total number of carbon atoms and the destruction rate of CO in molecules per second. Using the C mass range derived from our 1D model (section 4.1), we find a timescale between 1 500 and 8 000 years. When using the mass from the best fitting ‘ring + clump’ model of section 4.2.2, we find a timescale of 3 000 years. These are very short timescales compared to the age of the system. Thus, any event invoked to explain the observed CO clump (e.g. a giant collision) needs to occur at a correspondingly high rate—otherwise, it would be unlikely to observe the results of such an event so shortly after it occurred. In the following, we adopt 5 000 a as an average estimate of the production timescale.
A caveat remains that some of the C produced from CO photodissociation has already been removed. Kral et al. (2016) suggested that in steady state, atomic gas is removed by forming an accretion disk inside and a decretion disk outside of the CO-producing parent belt. However, our data argue against such a scenario (see section 4.3). Another possibility is chemical processes that might be able to change the amount of carbon in the disk (Higuchi et al., 2017). This would require a sufficiently high H density. Removing carbon by radiation pressure seems unlikely. First, both neutral and ionised C do not feel any radiation pressure around Pic. However, as shown by Fernández et al. (2006), ions are coupled into a single fluid via Coulomb interactions. But the effective radiation pressure on this fluid is not sufficient 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, accretion by a planet seems unlikely. Gap-opening planets (Jupiter-class) are the best candidates. But even such planets have a finite accretion efficiency (typically 75%–90%) limited by the leakage of flow across their gaps (Lubow & D’Angelo, 2006). So to increase the estimated lifetime by a factor of or more, the hypothetical planet would have to sit close to the observed belt and accrete with an efficiency of or higher. In addition, planets down to 1 M have been excluded 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 as estimated in section 5.3, the absence of an accretion disk is actually not surprising. Indeed, the timescale for viscous accretion from a radius can be estimated by
with the sound speed and the scale height. We assume with K and , and as in equation 2 with au. A very high (taking the upper bound of the timescale range calculated in section 5.3) would be required to match the viscous timescale with the time since C production started. In other words, for any reasonable value of , 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 Fe are seen in Keplerian rotation. These species are strongly affected by radiation pressure and C is needed to prevent them from being blown out. As was shown by Fernández et al. (2006), C, which does not feel any radiation pressure around Pic, is acting as a braking agent in the form of C via Coulomb interactions. In order to test whether the amount of carbon needed to brake metals is consistent with the new ALMA C i data, we consider an accretion disk model extending from 10 to 50 au with a carbon surface density and temperature (with K at 85 au) (Lynden-Bell & Pringle, 1974). We compute the emission from this model as described in section 4.2. We then search for C i emission only in those regions of the datacube where the model predicts emission. However, we want to exclude regions of the datacube that can contain emission from the outer disk. Thus, for those data points with (with the orbital speed at radius ), we request points to be at least one spatial resolution element inside of the line in the pv diagram defining a thin ring at 50 au (see figure 2). For points with a larger , we can be sure that no contamination from the outer disk is present. We also exclude points with a predicted emission below 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 what is described in appendix A. Defining a detection threshold of , we do not detect significant emission. We derive an upper limit on the C density by adjusting the model such that the probability for a non-detection, i.e. measuring a flux below the detection threshold, would only be 1% if the model was correct. Assuming Gaussian noise, the probability of a model with flux to remain undetected is given by with the error on the flux and 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ández et al., 2006). Thus, the data are consistent with enough carbon being present in the inner disk to brake metals. However, to get better constrains, we also look at the C ii observations from Herschel/HIFI by Cataldi et al. (2014). We measure the flux in the wings of the C ii spectrum corresponding to radial velocities between and . Errorbars are calculated by calculating the flux in spectral regions of the same size without line emission and taking the standard deviation. We first detemine the errors on the flux in the left 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 model such that the combined probability to remain undetected in the ALMA C i and the HIFI C ii (H and V beam) data is only 1% (including the ALMA data has a negligible effect as the HIFI data are more constraining). This model has a mid-plane C density at 30 au of 380 cm. We conclude that the currently available data are consistent with carbon being present inside of 50 au at a level that is sufficient to brake metals. We suggest that this gas in the inner region was 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 previous gas-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 the disk. This is surprising since one would expect C to spread in azimuth within a few orbits even though C production might happen preferably at the CO clump. This is in contrast to CO which has a lifetime shorter than an orbital period 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 started so recently that there was not yet enough time for azimuthal spreading respectively symmetrisation. We expect the timescale 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 the CO destruction rate calculated in section 5.3. In reality, only 30% of the CO emission is found in the clump (Dent et al., 2014). This setup thus maximises the asymmetry, so the derived symmetrisation timescale can be considered an upper limit. We then write the following simple equations describing the temporal evolution of the neutral and ionised carbon densities under the influence of C production, ionisation, recombination and orbital motion:
where and are the azimuthal densities (in kg rad) of neutral and ionised C respectively, is the CO destruction rate (in kg s), and are atomic and molecular masses respectively, the electron density, the recombination coefficient, the ionisation rate and the angular velocity at 85 au. For the electron density, we assume that all electrons come from the ionisation of C, and that the typical volume extends over 70 au in the radial direction (best-fit ‘ring’ model, section 4.2.1) and over (with the scale height, see equation 2) in the vertical direction (the model remains one-dimensional; we only assume a certain volume to get a value for the electron density). This means that where is the mass of a C ion. Then, equations 5 and 6 are used to compute the change of the azimuthal densities at each time step. Figure 13 shows how the SW/NE mass ratio of neutral carbon evolves with time. The wavy pattern is due to the orbital motion of the gas. The local minima correspond to the times when the first gas produced by the point source leaves the NE side and enters the SW side. After 10 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 we observe today if there is no mechanism that prevents symmetrisation. For example, assuming that only 30% of the C production occurs at the clump, with the other 70% uniformly distributed at all azimuths (this corresponds to adding a production term independent of to equation 5), we get the green curve shown in figure 13. As expected, the symmetrisation occurs faster—within 10 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 robust constraints, detailed hydrodynamical calculations would be necessary.
The upper limit on the symmetrisation timescale is below the lower bound of the timescale range over which C production occurred, as derived from the C mass (section 5.3). Thus, it appears unlikely that the asymmetry can be explained by invoking a very recent event.
In summary, the estimated lifetime (5 000 a) of the carbon-rich gas disk should be long enough to spread out the azimuthal asymmetry, but not long enough to diffuse the disk radially via viscous spreading.
5.5.2 Resonance trapping by planet
Matrà et al. (2017b) argued that the asymmetry in Pic, and perhaps in a number of other debris disks, is the result of planetesimals being trapped into mean-motion resonances (MMR) by a (migrating) planet. For Pic, our resolved C i observation excludes this possibility.
For planetesimals in resonance with a planet there would be an enhancement in particle density at some special azimuth relative to the planet. In the case of 2:1 MMR, such a resonance island lies behind and ahead the azimuth of the perturbing planet. As collisions are more frequent in the denser region, it is expected, in this scenario, that CO, which essentially traces recent collisions because of its short lifetime, is enhanced near the island and therefore asymmetric. However, as the planet revolves in its orbit, the resonance island sweeps through all azimuths. Integrated over many periods, the resonance island does not linger particularly long around a special azimuth. So if we look at a tracer gas that is the integrated production of CO over many periods, as the carbon gas is (produced over 5 000 a, i.e. 10 orbits), we expect it to be evenly spread out in the orbit. The observed asymmetry in C thus cannot be explained 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. One can show that the resonance width, even when forced by a highly eccentric Jovian planet, is narrower than the vertical scale height of the gas ( a few au). As the vertical scale height corresponds to the sound speed travel time over an orbit, gas pressure can, in an orbital period, easily disperse the gas azimuthally and radially over an extent wider than the 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 over a time much longer than the sound crossing time. However, over the lifetime of the system (20 Ma), differential precession should have disturbed this initial imprint. Ignoring the presence of other planets, and only considering the quadrupole precessional effect of Pic b (Lagrange et al., 2009), the timescale for order unity change in the relative precession angle is (Murray & Dermott, 2000)
where we evaluate the expression using M (Morzinski et al., 2015) for the planet mass and au (Wang et al., 2016) for its semi-major axis with the ring at au and with a fiducial radial width of au (consistent with that in Table 2). So even without an additional planet, an initially eccentric disk is expected to be largely 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. The age of the C gas (5 000 a) is shorter than any reasonable secular timescale (, where is the mass ratio of planet 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 a ago, if they have lived in such states 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 found near apoaps where they move the slowest (the so-called ‘apo-centre’ glow, Pan et al., 2016), while collisions are more likely to occur near the pericentre where the particle streamlines are more densely spaced and their relative velocities are higher. Collisions near the pericentre occur at a higher relative velocity, allowing smaller debris (which are more populous and have a larger collisional surface area) to break apart a given fragment. Matrà et al. (2017a) suggest that the mass loss rate is enhanced at either periaps or apoaps depending on the proper eccentricity and strength of the planetesimals. From preliminary numerical simulations, we favour an enhancement at periaps (more detailed simulations will be presented in a forthcoming paper). As a result, we expect CO (which reflects the instantaneous collisional mass loss rate) to be concentrated in periaps (the SW side, as is observed), while both the submm emission and the carbon line fluxes should be determined mostly by particle trajectories and should peak at apoaps, contrary to the observations (the observed C i periapsis to apoapsis flux ratio is , see section 3.1). This argument is easy to understand if one thinks of the CI gas as exact analogue of the small dust grains, which are also integrated products of previous collisions. Small grains will shine more brightly in apo-centre because there are more of them there—unless their eccentricities are so large that the fall-off of stellar-flux at apo-centre is more important than the number excess, see our discussion below.
At higher eccentricities, the apocentre is much further than the pericentre, and the drop-off in stellar flux, particle volume 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 the carbon gas is less excited, leading to reduced line fluxes. This effect is more severe when the orbits have an intrinsic spread in eccentricity, leading to a spread in apocentre distances and further reducing the dust and gas volume densities in those locations. Carbon ionization fraction, on the other hand, may also be affected by the lower ionizing flux and the smaller recombination rate.
So to explain the observed same-sided asymmetry in different tracers, we will require a medium to high eccentricity (, table 2) and preferably a large spread in eccentricities (discussed below). For a secularly forced ring, the forced eccentricity , where is the eccentricity of the planet (Murray & Dermott, 2000). A planet with considerable eccentricity is then required, which may itself decimate the planetesimal belt by dynamical ejection.
Another argument against the planet hypothesis comes from the range of eccentricities required to explain the observations. For a disk with a spread in semi-major axis ( au if au), we expect a spread in the forced eccentricities. This is too small to help explain the same-side asymmetry and it is also smaller than indicated by our best-fit solution.
In conclusion, it seems the hypothesis that an underlying highly-eccentric planetesimal belt, forced secularly by a planet, 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 produces a wide spread of debris that have a range of eccentricities but a similar alignment. This is qualitatively plausible to explain the observations (but see Matrà et al., 2017b, for a counter argument regarding the radial width). If such an event 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 the event rate.
Consider bodies with radius , in a belt of semi-major axis and width , performing vertical epi-cycles around the mid-plane. Viewed by each one of these bodies, the other bodies present a certain optical depth of , or a mean collision time of . Summed over bodies, this implies a mean event time of
where we assume there is no gravitational focusing among these low-mass objects, reasonable if their velocity dispersion lies much beyond their surface escape velocities. The scaled value km is appropriate to explain the total gas/dust mass observed, and AU is inspired by the best fit in Table 2. The value is a place-holder. So, to produce an event as recent as 5 000 a ago, one would require , or an absurdly high total mass of . This argument makes giant impacts exceedingly implausible at this location in the disk. Another difficulty with such a scenario is that giant impacts tend to be completely accretionary at low relative speeds, and even at speeds beyond the surface escape, only a small fraction of mass can be unbound and released into the circumstellar environment (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 be presented in an upcoming paper. Consider that the outer disk contains a number of Neptune-like planets ( a few). They have a surface escape velocity of km s. Let us also assume that there are bodies similar in mass to the Moon () or Mars () and moving in space with a relatively low dispersion velocity, . There is little gravitational focussing among these bodies, but strong focussing by the Neptunes. The timescale for a physical collision with a Neptune is
where we chose . This is somewhat arbitrary, but the abundance of cold Neptunes is already suggested by micro-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 a debris disk around the star. Instead, we focus on encounters that are close enough that the body is tidally disrupted. This means encounter distance of and the encounter time Ma (at , while the geometric cross-section goes up by a factor of , the gravitational focussing factor, , evaluated at goes down by a factor of ). So there could have been a few tidal disruption events over the lifetime of Pic, if there are thousands of moons floating around and if these moons retain low enough dispersion velocities to experience strong gravitational 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 will be ejected on a variety of orbits, bringing about a large spread in eccentricity. The semi-major axis will also have a spread. All debris will return to the disruption site, making this location a region of frequent collisions. The size of the disruption site, however, is larger than that in giant impact. As the debris flies away from Neptune, its stellar-centric orbit is altered by Neptune’s gravity while it is still within Neptune’s Hill sphere. As a result, unlike the narrow nozzle of the size of the impactor radius in the case of giant impacts, here, the nozzle has a typical spread of order the planet’s Hill radius, which reduces the peak collision rate. However, Matrà et al. (2017b) measured a radial extent of the CO clump of at least 100 au, which is clearly larger than the planet’s Hill sphere. More detailed modelling is needed to see whether a tidal disruption event can produce such a radially broad clump. Also, the clump in the deprojection of Matrà et al. (2017b) might appear more extended than it is in reality because the intrinsic velocity dispersion of the gas and the finite velocity resolution of the instrument degrade the resolution of the deprojection in the direction. In addition, the deprojection is carried out assuming velocities corresponding to circular orbits. Thus, if the orbits are eccentric 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 special time. The duration of a tidal disruption flare will depend on the above-discussed collisional geometry, but also on the distribution of the largest fragments (that remain or reform) after the disruption event. Both of these effects need to be studied in detail. In addition, a short lifetime will also ensure that debris products from previous events do not interfere with the current event. If not, the asymmetry would be washed out by previous debris which likely have a different 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, secular forcing). Instead, we propose that the bright, asymmetric debris disk in Pic could be the result of a recent tidal disruption of a Moon to Mars-sized object by a Neptune-like planet.
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 spatial and spectral distribution of the emission really is. To answer this question, we first interpolate the CO data cube onto the coordinates of the C i data cube. Then, we use convolution to adjust the spatio-spectral resolution of the data cubes. Finally, we normalise the data cubes and subtract. Figure 14 shows the moment 0 and pv diagram of the residual cube for the CO 2–1 and 3–2 transitions (Matrà 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 in the emission distribution. Note that there is significant difference in the distribution of CO 2–1 vs CO 3–2 emission (Matrà 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 of what we see for C and CO (compare figure 3 of this work to figure 3 in Brandeker et al., 2004). Na and Fe also show a tilt similar to what is seen for C and CO. The radial distribution is quite different, however, with the density being higher closer to the star instead of peaking at 85 au. The asymmetry is seen in both the inner and outer parts of the disk; concerning the inner distribution of Na and Fe, the emission is traced inwards to the observational limit of 13 au in the NE, while to the SW the density peak appears to be located much further out, beyond 50 au. In the outer parts of the disk, Na and Fe can be seen all the way to the limit of the observation at the projected distance of 330 au in the 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 the origin of the metals observed in the optical (Na, Fe, Ca, Ni, Ti, and Cr, Brandeker et al., 2004). While CO likely comes from the disruption of CO-rich cometary bodies at 85 au, the metals observed in the optical could be produced by the so-called falling evaporating bodies (FEBs, e.g. Vidal-Madjar et al., 2017) close to the star and then diffused outwards by the stellar radiation pressure. As the presence of C would act as a braking agent (Fernández et al., 2006), its spatial distribution would naturally become imprinted on the spatial distribution of the metals. The spatial distribution of Na 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 of eccentric orbits, the C would be more concentrated to the clump location in the SW, while being much more radially distributed in the NE, in agreement with the observed Na and Fe distributions.
6 Summary and conclusions
In this paper, we present resolved ALMA band 8 observations of C i emission towards the Pic debris disk. Our work can be summarised as follows:
Using a simple 1D model that calculates the ionisation balance and non-LTE level populations, we estimate the total C mass to be between and M. Assuming that C is produced from the photodissociation of CO at a constant rate, and that C is not removed from the system, this mass implies that gas production started only 5 000 a ago.
Surprisingly, C i shows the same asymmetry as seen for CO: a clump on the SW side of the disk. By modelling the spatial distribution of the C gas, we find that a satisfactory fit to the C i and archival C ii data can be found by assuming that the gas consists of a ring between 50 and 120 au with a superimposed clump at the same location as the CO clump. A model assuming eccentric orbits of the gas with a flat eccentricity distribution between 0 and 0.3 also reasonably fits the data.
The C i data are not consistent with the accretion disk predicted by Kral et al. (2016). If C gas production indeed started only 5 000 a ago, not enough time has passed for gas to have spread viscously into an accretion disk.
The short timescale since gas production started argues against a giant impact-origin of the C/CO/dust clump, because giant impacts are rare. It is unlikely that we observe the results of such an event so shortly after it occurred. However, while the production timescale of 5 000 a is short compared to the age of the system, it is long enough to allow azimuthal spreading of the gas. Thus, a scenario where the C assymetry is due to a lack of time for azimuthal spreading is disfavoured.
The fact that C shows the same asymmetry as CO (a clump in the SW) argues against a scenario where the clump is due to planetesimals trapped in a resonance with an outward-migrating planet. Indeed, in such a scenario, the clump orbits with the planet, i.e. the gas production should be symmetric on an orbital timescale.
In order to explain the simultaneous C and CO asymmetry, we propose that the planetesimal and gas disk of Pic is eccentric and might have originated from a recent tidal disruption event. This could potentially also explain 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 observations of the C i line at higher SNR and spectro-spatial resolution can be used to confirm or reject our hypothesis of an eccentric disk due to a tidal disruption event.
Appendix A Details of the error calculation
The total line emission was measured by integrating the data cube within 6 km s in the spectral dimesion and over the rectangular box shown in Figure 1 in the spatial dimension. Having measured the noise in the data cube, a naive way to estimate the error on the total flux would be to write with the number of data points over which the integration extends, the noise in the data cube and and the extent in solid angle and velocity 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 an estimate of the error. The volumes are placed such that they are sufficiently distant from each other to be considered independent. Since no primary beam correction has been applied, the noise can be assumed uniform over the data cube.
For the continuum, the above procedure yields too few flux samples. Thus, we reduce the size of the flux samples in the direction to get more samples. Then, we set where is the standard deviation of the flux samples and is the number of flux samples that fit into the region for which the flux is measured. Again, the flux samples are 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 one pixel in the direction. In the case of the profile derived from the restricted region in pv space (Figure 3), the size of the region over which we integrate depends on . Thus, we take flux samples for each individually. For this profile, the error thus depends on .
Appendix B Measurement of the line broadening parameter
The line broadening parameter (defined as where is the full width at half maximum of the line) parametrises the line broadening due to effects other than the orbital motion of the gas (e.g. thermal broadening or turbulence). We can use the ALMA observations of C and CO to measure by considering a vertical cut in the PV diagram (Figure 2) going through , i.e. considering the line of sight towards the star. For circular orbits, all gas along this line of sight is centred at the same radial velocity (namely 0 km s), allowing a measurement of . However, if orbits are eccentric, the line of sight towards the star can contain additional broadening due to orbital motion. Furthermore, the finite spatial resolution of the instrument will blur material with non-zero radial velocity into the line of sight towards the star. Thus, our derived value of should be considered an upper limit.
We consider the PV diagrams of C i (this work), CO 2–1 and CO 3–2 (Matrà et al., 2017b) and for each of them fit a Gaussian to the vertical cut through , yielding three independent measurements of where is the broadening parameter of the fitted Gaussian and is the broadening due to the spectral response function of the instrument777See ALMA Cycle 6 Technical Handbook, section 5.5.2 (Spectral Setup), https://almascience.nao.ac.jp/documents-and-tools/cycle6/alma-technical-handbook. We estimate errors using emcee, a python implementation of an affine invariant Markov chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al., 2013). The likelihood was defined via a measure analogous to equations 1 and 2 by Booth et al. (2017). In particular, we also use a noise correlation ratio, defined in our case as the square root of the number of spectral pixels per spectral resolution element. We consider flat priors, imposing only that the peak of the Gaussian and the broadening parameter are positive. We use 100 walkers with steps. This yields the following results for : km s (C i), km s (CO 2–1) and km s (CO 3–2). Combining these measurements as a weighted mean yields km s. This value may be compared to previous measurements of . Jolly et al. (1998) measured km s from CO absorption (they also measured a high value of 4.2 km s from C i absorption, potentially caused by modelling difficulties because the line is saturated). Roberge et al. (2000) obtained km s from C i absorption and km s 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 value derived from the ALMA data seems thus generally lower than in previous publications.
Appendix C Details of the deprojection procedure
The pv diagram can be used to get a deprojected view of the emission in the plane. As is discussed in appendix C of Matrà et al. (2017b), a radius can be assigned to points of the pv diagram for an edge-on disk and circular Keplerian orbits. From this, we can find the position along the line of sight . Note that for some points of the pv diagram, the term under the square root becomes negative. This simply means that for the given , the radial velocity cannot be reached anywhere along the line of sight, i.e. we have , where is the orbital velocity of the orbit with . One has to choose how to distribute the flux from a given pv point among the two possible 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 ). Then, for a given point , we calculate . From this, , where the minus sign accounts for the known rotation sense of the Pic disk. We then look up the flux at in the pv diagram and assign it to the point in the deprojection.
In order to get a deprojection with a straightforward interpretation in terms of the distribution of the flux, it is also 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 to . The Jacobian determinant of this transformation (for as in figure 4) reads
The flux read from the pv diagram is then multiplied by and an additional constant factor (with the central frequency of the line and the distance between the observer and Pic).
If the flux units of the deprojection are not transformed (as for example in Matrà et al., 2017b), an interpretation of the deprojection is not straightforward because the flux in a certain area of the deprojection does not equal the integral over and over this area. Figure 15 shows the deprojection without multiplication by the Jacobian (i.e. in the same units as the pv diagram). Comparing to figure 4, it becomes apparent that a deprojection without transformed units might lead to visual mis-interpretation, for example about the relative amount of flux at large radii (say beyond 150 au). Indeed, from equation C1, we see that the Jacobian gets smaller at large radii.
Appendix D Approximate calculation of the statistical equilibrium
We 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 dust continuum) is not attenuated by the gas. In this section, we justify these assumptions.
The SE equation for an atomic level can be written
where is the fractional population of level , is the collisional (de)excitation rate (with the collisional rate coefficient and the electron density), the frequency-integrated mean intensity and and are the Einstein coefficients (where we set if ). If we can show that (for any ), then we have 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 of the disk: . To get insight into the individual contribution of the two components, we consider them separately. First, we calculate for all transitions and all locations (except where the gas 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. We find that for C i and for C ii. Thus, although we included background radiation in our calculation, it is actually negligible. We double-check by re-computing the models without background radiation and find that the results indeed do not change.
Next, we consider (de)excitation by line emission and compute . 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 the other grid points, neglecting optical depth and the velocity field (the velocity field could red/blue-shift photons from other grid points out of the transition). We are thus calculating an upper limit on . We find that for C i and for C ii. Thus, a more sophisticated model should include the full radiative transfer for C i, but neglecting the line photons is a decent approximation given the quality of our data, since it considerably simplifies the calculation of the models.
As an additional test, we used the LIME code version 1.9.1 (Brinch & Hogerheijde, 2010) to calculate the full non-LTE radiative 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.25 (for C i) respectively 1.18 (for C ii).
Appendix E Gas density for eccentric orbits
In this section, we derive the gas surface density of the model with eccentric orbits presented in section 4.2.3. We assume that all orbits have a common pericentre and that the distribution of eccentricities is known and given by . We search the probability (which we assume is proportional to the surface density) to find a particle at radius and true anomaly . We first consider . In our model, a given eccentricity corresponds to a single orbit (because all orbits have a common pericentre). Thus, is interpreted as the probability to find a particle at true anomaly along the orbit with eccentricity . The time a particle spends at a given is inversely proportional to the orbital velocity. Thus, we have
where the normalisation constant . The orbital velocity is given by
with and the fixed pericentre distance. Next, we consider the transformation from to where
The transformation is bijective, except for . Therefore, we can write
with the Jacobian determinant of the transformation, given by , which can be calculated from equation E3. At , the density diverges. In practice, this is easily handled by simply not sampling the singularity on the numerical grid.
- Absil et al. (2013) Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12
- Agnor & Asphaug (2004) Agnor, C., & Asphaug, E. 2004, ApJ, 613, L157
- Apai et al. (2015) Apai, D., Schneider, G., Grady, C. A., et al. 2015, ApJ, 800, 136
- Armitage (2009) Armitage, P. J. 2009, Astrophysics of Planet Formation (Cambridge University Press), doi:10.1017/CBO9780511802225
- Artymowicz (1997) Artymowicz, P. 1997, Annual Review of Earth and Planetary Sciences, 25, 175
- Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
- Aumann (1985) Aumann, H. H. 1985, PASP, 97, 885
- Backman & Paresce (1993) Backman, D. E., & Paresce, F. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine, 1253–1304
- Booth et al. (2017) Booth, M., Dent, W. R. F., Jordán, A., et al. 2017, MNRAS, 469, 3200
- Bouret et al. (2002) Bouret, J.-C., Deleuil, M., Lanz, T., et al. 2002, A&A, 390, 1049
- Brandeker (2011) Brandeker, A. 2011, ApJ, 729, 122
- Brandeker et al. (2004) Brandeker, A., Liseau, R., Olofsson, G., & Fridlund, M. 2004, A&A, 413, 681
- Brandeker et al. (2016) Brandeker, A., Cataldi, G., Olofsson, G., et al. 2016, A&A, 591, A27
- Brinch & Hogerheijde (2010) Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25
- Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167
- Cataldi et al. (2014) Cataldi, G., Brandeker, A., Olofsson, G., et al. 2014, A&A, 563, A66
- Chilcote et al. (2017) Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182
- Crawford et al. (1994) Crawford, I. A., Spyromilio, J., Barlow, M. J., Diego, F., & Lagrange, A. M. 1994, MNRAS, 266, L65
- Crifo et al. (1997) Crifo, F., Vidal-Madjar, A., Lallement, R., Ferlet, R., & Gerbaldi, M. 1997, A&A, 320, L29
- Dent et al. (2014) Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Science, 343, 1490
- Draine (1978) Draine, B. T. 1978, ApJS, 36, 595
- Fernández et al. (2006) Fernández, R., Brandeker, A., & Wu, Y. 2006, ApJ, 643, 509
- Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
- Gray et al. (2006) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161
- Grigorieva et al. (2007) Grigorieva, A., Thébault, P., Artymowicz, P., & Brandeker, A. 2007, A&A, 475, 755
- Heays et al. (2017) Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105
- Higuchi et al. (2017) Higuchi, A. E., Sato, A., Tsukagoshi, T., et al. 2017, ApJ, 839, L14
- Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
- Jackson et al. (2014) Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D. 2014, MNRAS, 440, 3757
- Jolly et al. (1998) Jolly, A., McPhate, J. B., Lecavelier, A., et al. 1998, A&A, 329, 1028
- Kiefer et al. (2014) Kiefer, F., Lecavelier des Etangs, A., Boissier, J., et al. 2014, Nature, 514, 462
- Kral & Latter (2016) Kral, Q., & Latter, H. 2016, MNRAS, 461, 1614
- Kral et al. (2017) Kral, Q., Matrà, L., Wyatt, M. C., & Kennedy, G. M. 2017, MNRAS, 469, 521
- Kral et al. (2016) Kral, Q., Wyatt, M., Carswell, R. F., et al. 2016, MNRAS, 461, 845
- Lagrange et al. (2009) Lagrange, A.-M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21
- Lagrange et al. (2010) Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57
- Lagrange et al. (2012) Lagrange, A.-M., Boccaletti, A., Milli, J., et al. 2012, A&A, 542, A40
- Lecavelier des Etangs & Vidal-Madjar (2016) Lecavelier des Etangs, A., & Vidal-Madjar, A. 2016, A&A, 588, A60
- Lee (1984) Lee, L. C. 1984, ApJ, 282, 172
- Li et al. (2012) Li, D., Telesco, C. M., & Wright, C. M. 2012, ApJ, 759, 81
- Lieman-Sifry et al. (2016) Lieman-Sifry, J., Hughes, A. M., Carpenter, J. M., et al. 2016, ApJ, 828, 25
- Liseau (2003) Liseau, R. 2003, in ESA Special Publication, Vol. 539, Earths: DARWIN/TPF and the Search for Extrasolar Terrestrial Planets, ed. M. Fridlund, T. Henning, & H. Lacoste, 135–142
- Lubow & D’Angelo (2006) Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526
- Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
- Mamajek & Bell (2014) Mamajek, E. E., & Bell, C. P. M. 2014, MNRAS, 445, 2169
- Matrà et al. (2018) Matrà, L., Wilner, D. J., Öberg, K. I., et al. 2018, ApJ, 853, 147
- Matrà et al. (2017a) Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017a, ApJ, 842, 9
- Matrà et al. (2017b) Matrà, L., Dent, W. R. F., Wyatt, M. C., et al. 2017b, MNRAS, 464, 1415
- Matthews et al. (2014) Matthews, B. C., Krivov, A. V., Wyatt, M. C., Bryden, G., & Eiroa, C. 2014, Protostars and Planets VI, 521
- McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
- Millar-Blanchaer et al. (2015) Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, 811, 18
- Milli et al. (2014) Milli, J., Lagrange, A.-M., Mawet, D., et al. 2014, A&A, 566, A91
- Moór et al. (2011) Moór, A., Ábrahám, P., Juhász, A., et al. 2011, ApJ, 740, L7
- Morzinski et al. (2015) Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108
- Mumma & Charnley (2011) Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471
- Murray & Dermott (2000) Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics
- Nesvold & Kuchner (2015) Nesvold, E. R., & Kuchner, M. J. 2015, ApJ, 815, 61
- Nilsson et al. (2012) Nilsson, R., Brandeker, A., Olofsson, G., et al. 2012, A&A, 544, A134
- Olofsson et al. (2001) Olofsson, G., Liseau, R., & Brandeker, A. 2001, ApJ, 563, L77
- Pan et al. (2016) Pan, M., Nesvold, E. R., & Kuchner, M. J. 2016, ApJ, 832, 81
- Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
- Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2
- Redfield (2007) Redfield, S. 2007, ApJ, 656, L97
- Roberge et al. (2000) Roberge, A., Feldman, P. D., Lagrange, A. M., et al. 2000, ApJ, 538, 904
- Roberge et al. (2006) Roberge, A., Feldman, P. D., Weinberger, A. J., Deleuil, M., & Bouret, J.-C. 2006, Nature, 441, 724
- Roberge et al. (2014) Roberge, A., Welsh, B. Y., Kamp, I., Weinberger, A. J., & Grady, C. A. 2014, ApJ, 796, L11
- Rollins & Rawlings (2012) Rollins, R. P., & Rawlings, J. M. C. 2012, MNRAS, 427, 2328
- Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
- Slettebak (1975) Slettebak, A. 1975, ApJ, 197, 137
- Slettebak & Carpenter (1983) Slettebak, A., & Carpenter, K. G. 1983, ApJS, 53, 869
- Smith & Terrile (1984) Smith, B. A., & Terrile, R. J. 1984, Science, 226, 1421
- Telesco et al. (2005) Telesco, C. M., Fisher, R. S., Wyatt, M. C., et al. 2005, Nature, 433, 133
- van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
- van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
- Vandenbussche et al. (2010) Vandenbussche, B., Sibthorpe, B., Acke, B., et al. 2010, A&A, 518, L133
- Vidal-Madjar et al. (2017) Vidal-Madjar, A., Kiefer, F., Lecavelier des Etangs, A., et al. 2017, ArXiv e-prints, arXiv:1709.08170
- Vidal-Madjar et al. (1994) Vidal-Madjar, A., Lagrange-Henri, A.-M., Feldman, P. D., et al. 1994, A&A, 290, 245
- Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323
- Wang et al. (2016) Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97
- Welsh & Montgomery (2015) Welsh, B. Y., & Montgomery, S. L. 2015, Advances in Astronomy, 2015, 980323
- Wilson et al. (2017) Wilson, P. A., Lecavelier des Etangs, A., Vidal-Madjar, A., et al. 2017, A&A, 599, A75
- Wyatt (2003) Wyatt, M. C. 2003, ApJ, 598, 1321
- Wyatt (2006) —. 2006, ApJ, 639, 1153
- Wyatt (2008) —. 2008, ARA&A, 46, 339
- Xie et al. (2013) Xie, J.-W., Brandeker, A., & Wu, Y. 2013, ApJ, 762, 114
- Zagorovsky et al. (2010) Zagorovsky, K., Brandeker, A., & Wu, Y. 2010, ApJ, 720, 923