Exocomets size distribution in the [Formula: see text] Pictoris planetary system.

The star [Formula: see text] Pictoris harbors a young planetary system of about 20 million years old, which is characterized by the presence of a gaseous and dusty debris disk, at least two massive planets and many minor bodies. For more than thirty years, exocomets transiting the star have been detected using spectroscopy, probing the gaseous part of the cometary comas and tails. The detection of the dusty component of the tails can be performed through photometric observations of the transits. Since 2018, the Transiting Exoplanet Survey Satellite has observed [Formula: see text] Pic for a total of 156 days. Here we report an analysis of the TESS photometric data set with the identification of a total of 30 transits of exocomets. Our statistical analysis shows that the number of transiting exocomet events (N) as a function of the absorption depth (AD) in the light curve follows a power law in the form [Formula: see text], where [Formula: see text]. This distribution of absorption depth leads to a differential comet size distribution proportional to [Formula: see text], where [Formula: see text], showing a striking similarity to the size distribution of comets in the Solar system and the distribution of a collisionally relaxed population ([Formula: see text]).

Introduction

Since the mid-1980’s, spectroscopic observations of the bright southern star  Pictoris, or  Pic, have revealed variations in the calcium H and K lines, which have been interpreted as the transit of the gaseous tails of exocomets[1-5]. Several decades of exocomets transit observations, mostly on  Pic but also on a few other stars, provided original insight on the activity of minor bodies in planetary system within the first tens of million years[6]. For that, the  Pic system is exceptional in presenting simultaneously several favorable properties : it is only a few tens million years old, very nearby at 19.3 parsecs and seen exactly edge-on from the Earth[7-16]. Nonetheless, even for  Pic, the observations of exocomets were limited to the gaseous part of the tails probed by transit spectroscopy; observations of the dusty part of the tails indeed require high accuracy photometric measurements, which was until recently not available[17-21].

The situation changed when, using data obtained by the Transiting Exoplanet Survey Satellite (TESS[21]) between October 2018 and February 2019, three photometric events were discovered[20] and interpreted as due to the transit of the dust component of exocomets transiting  Pic. In support of this interpretation, the observed light curves are almost identical to the predictions made twenty years earlier[17,18].

Since 2019,  Pic has been re-observed by TESS. Here we present the analysis of the complete data set gathered up to February 2021 in order to perform a deep search for exocomet transits and determine the size distribution of the  Pic comets to be compared with the distribution observed in the Solar system.

Search for exocomet transits

Considering all TESS observations of  Pic until February 2021, after the cleaning process of the Scuti variations and other slower variations (see Methods) we end up with a 156 days light curve clearly showing some dips, which are similar to what is expected for exocomet transits (Extended Data Fig. 4). To make sure that the observed dips in the  Pic light curve are real and not noise residuals nor artifacts due to the reduction process, we used the TESS observations of the nearby star  Pic. We downloaded the PDC flux time series of  Pic and applied the same procedure to remove the Scuti and slower variations. The observations of  Pic provide an excellent data-set to test our procedure because  Pic is in the same region of the sky as  Pic (and has hence overlapping TESS observations epochs), it has the same spectral type (A8V versus A5V for  Pic), similar magnitude (3.3 versus 3.85 for  Pic) and a similarly high projected rotational velocity (205 km/s for  Pic versus 139 km/s for  Pic). Thus,  Pic is almost a nearby stellar twin of  Pic except for the presence of the young planetary system.  Pic has already been successfully used as a reference star for analysis of  Pic observations[22].

The  Pic light curve is used to check for systematics that could mimic transit of exocomets. The light curve of  Pic shows noisy excursions from the mean value both in the positive and negative deviations with the same pattern. The light curve of  Pic also shows noisy excursions in the positive direction that are similar to the ones observed in the light curve of  Pic, but noisy excursions that are more pronounced and more frequent in the negative direction. The latest are typical signatures of small, transiting objects with extended dust tails. We do not detect similar variations in  Pic that the ones seen in  Pic that can be due to the passage of exocomets in front of the star.

To characterize this excess of transit-like features in the light curve of  Pic and identify the corresponding individual events, we calculated the correlation of the light curves with a simple model of an exocomet transit photometric event, assuming a 1D transit of a translucent dust cloud with an exponential decrease of the optical thickness from the head of the comet. This model has only four parameters : K, the cloud optical thickness at the leading head, , the transit duration corresponding to the time needed to cover the chord length of the stellar disk at the transit velocity, , the speed of the transit of one scale length of the cometary tail (the inverse of the scale length of the cometary tail divided by the transit velocity), and finally , the time of the beginning of the transit. With this model of the transit of the exocometary tail, the corresponding decrease in relative flux at the time t is given by , where if , if , and if , and if .

We calculated the correlation of this 1D-model with the observed light curves of  Pic and  Pic by varying the value of the transit time and using various plausible values for the and parameters characterizing the shapes of the exocomet transit light curves. We used from 2 to 20  and from 0.15 to 0.5 days, corresponding to periastron distances ranging from about 0.08 to 0.85 au. The value of K in the model can be arbitrarily chosen because it only changes the amplitude of the light variations, hence it has no consequence in the position of the correlation maximum nor on the identification of exocomet transits events.

In the case of  Pic the correlation reaches large positive values for some of the values of . This behaviour is not observed with the  Pic light curve, showing that the photometric variations with cometary transit shapes are specific to  Pic. We check the negative values of the correlation of the model with the  Pic light curve itself, and found that they are much less numerous and significantly smaller than the positive values. This confirms that the light curve of  Pic shows photometric variations with decrease of the star brightness that are typical of the transits of exocomets and the absence of variations with increase of the star brightness with similar shape and amplitude. We interpret the correlation peaks as the signatures of potential exocomet transits. We keep only the peaks with correlation values that are higher than the maximum value obtained with the  Pic data and the maximum negative value with the  Pic data. This is our conservative criterion to consider the observed variations as a detection of an exocomet transit.

We identified a total of 30 significant detections of exocomet transit events. Although this is not surprising, acknowledging the well known ubiquity of comets in that young planetary system, this is the first time such a large number of exocomets are detected in photometry. This allows statistical analysis of their properties. The light curves of these detected exocomet transits are plotted in Extended Data Fig. 4. Their characteristics and the parameters of the best fits with the 1-D model are given in Extended Data Table 2. The last column of the table gives the square root of the improvements of the fits to the light curve between a model with no transit and the 1-D transit model using the best fit parameters. This shows that the identified exocomet transits are all detected at least at 4- level.

Size distribution of  Pictoris exocomets

Distribution of absorption depths

Our detection of 30 photometric transits of exocomets allows a statistical analysis of their properties. Here we call “absorption depth”, noted AD, the decrease in relative flux at the minimum of a transit light curve. The numerical value of the absorption depth is estimated from the best fit with the 1-D model with . A plot of the events frequency as a function of the absorption depth shows that there is a steep decrease of the number of events toward the larger absorption depths (Fig. 1). The differential number of transiting exocomet events (dN) detected with an observation of duration as a function of the absorption depth can be fitted by a power law in the form . Considering 29 events detected in 156 days of observations, that is all the 30 events except the deepest one (see below), we find that and , where the uncertainties have been evaluated using a Poisson distribution for the number of events in each bin of width .

Figure 1

Histogram of the number of exocomet transit events as a function of the absorption depth. The red line shows the fit with a power law function with . The uncertainty on the fitted parameters of the power law function have been evaluated using a Poisson distribution for the number of events in each bin of width . The red squares represents the number of expected events in each bin as calculated with the fitted power law.

The deepest transit event[20] of Julian Day JD = 2457000 + 1486 looks exceptional with an absorption depth of about . This event could be produced by a member of another family of exocomets than the one which produces the 29 other shallower events, as we know from spectroscopic transit observations the presence of several families of  Pic comets[5]. Nonetheless, with the distribution derived above, the expected number of events with absorption depths in the range [10-20] in a 156 days observation is . Even in the range [15-20] the expected number of detections in 156 days is , corresponding to a probability of 26% to have one single event in this range as observed. Therefore the event of JD = 2457000 + 1486 can simply be a normal event, only the deepest, within the same distribution of the other 29 events. If this event is taken into account, the estimates of the distribution power law indexes and change by about 0.1 and 0.2, towards slightly shallower distributions. Nonetheless, to remain conservative, this deepest event is not taken into account in the derivation of the distribution of absorption depths and exocomets sizes considered here.

Distribution of exocomet sizes

The modeling of the exocomet transit light curves shows that the transit absorption depth, AD, is directly proportional to , the dust evaporation rate from the comet nucleus[17] . If we assume that the dust production rate is proportional to the comet nucleus area (i.e., Ref.[23]), we have a production rate proportional to the squared of the nucleus radius R. Finally, with AD proportional to , we find that the differential number of exocomets as a function of the nucleus size is given by , with .

With the fit to the observed distribution of the absorption depths, we conclude that the differential distribution of the exocomet size must follow a power law with an index . This distribution is notably similar to the size distribution of comets in the Solar system (Fig. 2) and the distribution predicted in Ref.[24] for a collisionally relaxed population ().

Figure 2

Plot of the cumulative size distribution of the exocomets in  Pic. The cumulative size distribution is plotted with blue squares for each exocomet and the corresponding fit excluding the largest comet is plotted with the red thick line. For comparison, published size distributions measured in the Solar system are plotted with thin dashed lines for asteroids in comets orbits (ACO), on near Earth orbits (NEO) and non-near Earth orbits (non-NEO) (A06, Ref.[31]), Jupiter family comets (JFC) (T06, Ref.[26]; S11, Ref.[27]; F13, Ref.[29]; B17, Ref.[30]), and long-period comets (LPC) (B17, Ref.[30]; B19, Ref.[25]). In this plot, the size distributions for the objects in the Solar system have been scaled to have a cumulative number of about 10 objects with radius above 2 kilometers. The radii of the  Pic comets have been estimated using the derivation described in the Method section. The conclusion on the similarity of the size distributions in  Pic and the Solar system is independent of these estimates.

For the plot of the size distribution (Fig. 2), we used the cometary radii estimated following the derivation described in the Method section. The conclusion on the similarity of the size distributions in  Pic and the Solar system is independent of these absolute size estimates. Nonetheless, it is remarkable that not only the distribution but also the sizes of the  Pic comets nuclei are found to be similar to sizes of the Solar system comets.

Comparison with Solar system comets

The size distribution of the nucleus of Solar system comets has been estimated for various locations. The size distribution of the Jupiter family comets and Oort cloud comets are found to be similar but not exactly the same.

For the Oort cloud, the size distribution has been estimated using a cometary activity model with a survey simulation and application to 150 long-period comets (LPC) detected over 7 years by the Pan-STARRS1 near-Earth object survey[25]. For objects with diameters above 1 km, the distribution is found to be  (stat.)  (sys.), which is in remarkable agreement with our value for  Pic exocomets. For smaller long-period comets, a shallower distribution is found with for diameters between 100 m and 1 km (Ref.[25]). Note that in practice, in a similar manner as we have used the transit absorption depth as a proxy for the size estimate of the  Pic comets, the above estimates for the Oort cloud comets have been obtained using the absolute H magnitudes of the nuclei as the proxy for their size.

For the  Pic exocomets, the comparison may be more appropriate with the size distribution of the comets in the Jupiter family, which originates from the Kuiper Belt. In Ref.[26] a catalog of absolute nuclear magnitudes of Jupiter family comets (JFC) has been used to derive a size distribution and to find for nuclei with radius between 2 and 5.5 kilometers (see Fig. 8 of Ref.[26]). More recent works have provided a shallower distribution : analyzing a large number of optical observations, Ref.[27] found a lower value with for nuclei with radius larger than 1.25 kilometers, and for smaller objects. This last result is consistent with the result described in Ref.[28], where images of Jupiter family comets obtained with the Hubble space telescope and the Keck telescopes have been analyzed. With a model fit to the observations, it is concluded that the intrinsic size distribution of comets in the Jupiter family is consistent with a power-law but truncated at small nucleus radii below 2.0 kilometers. In Ref.[29] a similar distribution is obtained with for radius between 2 and 5 kilometers, interpreting the distribution shallower than the canonical Dohnanyi’s size distribution (Ref.[24]) as due to fragmentation of the JFC objects. In Ref.[30] the measured for Jupiter family comets of radius between 2 and 10 kilometers is to be compared to the for long-period comets between 1 and 20 kilometers in radius.

Finally, in the Solar system the size distribution of extinct or dormant comets can be determined through the population of asteroids in comet orbits (ACO). Assuming that the exhaustion processes have no major impact on the size distribution, the observed distribution for these asteroids can be considered as a good proxy for the distribution for the parent comets. For these objects, Ref.[31] found for the full sample with radius between 2.8 and 7 kilometers, for near Earth objects (NEO) with radius down to 1.4 kilometers and for non-near Earth objects (non-NEO) with radius between 2.8 and 7 kilometers.

Taken all together these estimates for the Solar system comets are in general agreement with the value that we obtained for the  Pic comets, with some slightly shallower distributions in some cases (Fig. 2). This points toward the importance of collisional fragmentation in shaping the size distribution of the exocomets in the younger  Pic planetary system.

Discussion

The measured absorption depth distribution is the result of the distribution of several parameters for each individual comet, e.g., the orbital parameters, cometary activity, composition, size, etc. Here, following the result of numerical simulations[17,18], we assumed that the absorption depth distribution is mainly dominated by the distribution of the exocomet intrinsic dust production rate and hence their size. In other terms, although other parameters play a role for each individual comet, their diversities are expected to have a lower impact on the observed transit absorption depth than the size. In support of that idea, in spectroscopy it is observed that the  Pic transiting comets present similar orbital characteristics, which allows the classification in two different families[5]. With similar orbits, different transiting exocomets have different dust tails mainly because of different dust production rate, and hence because of different size nuclei.

The observed distribution of exocomets in the young planetary system of  Pic is strikingly similar to the distribution observed in the Solar system. This distribution seems to be ubiquitous and is also consistent with the canonical Dohnanyi’s size distribution[24] (), which corresponds to the size distribution of a collisionally relaxed population (see discussion in Ref.[32]). This indicates that the collisional process with fragmentation cascades is likely one of the dominant processes that shape the population of kilometer-sized bodies in the  Pic planetary systems.

Methods

TESS observations

Pictoris has been observed by TESS at 2-minute cadence several times from October 2018 to February 2021. The available data-set covers a total of 156 days of observations in the optical domain, divided into 14 epochs of about 12 days each. The raw data shows a flux dispersion of about , which is mainly due to Scuti pulsations in the stellar atmosphere.

The data from 19 October 2018 to 1 February 2019 have already been analysed[20]. Three photometric events have been identified and attributed to the transits of three different exocomets, with one spectacular transit at Julian Day (JD) equals to (2457000+1486).

In addition to these pioneering observations,  Pic has been observed from 20 November 2020 to 8 February 2021. For the observations of Sectors 4 to 7 and Sectors 31 to 34 (Extended Data Table 1), the TESS data products were obtained from the Mikulski Archive for Space Telescopes (MAST). We used the Presearch Data Conditioning (PDC) flux time series[33-35] processed by the TESS Science Processing Operations Center (SPOC) pipeline version 5.0.23-20210212[36,37], with the spacecraft cosmic ray mitigation algorithm enabled.

Cleaning the light curves

The Scuti variations

Pictoris is prone to Scuti type photometric variations[38-40]. These variations have a dominant frequency of (corresponding to a period of about 0.5 hours), with an amplitude of up to (Extended Data Fig. 1). These variations are superimposed to the photometric signatures of exocomet transits and must be corrected before searching for these exocomets.

Each of the 14 epochs of continuous observations (Extended Data Table 1) have been reduced separately; for each of them we extracted the set of frequencies and amplitudes of the pulsations using the Period04 software as previously done in Ref.[20]. The Period04 software performs a Fourier transform on the data, giving the frequencies, amplitude, phase and signal to noise ratio of each harmonic in the time series. We conducted the frequency search between 15 and , to avoid removing slow variations that might be caused by a cometary transit. For each iteration of the software, the highest-amplitude frequency within the search range is selected, and added to a multi-sine model. This model is then optimized over amplitude and frequency of each harmonic, then removed from the original signal. For each iteration, the software computes the signal to noise ratio of the main remaining frequencies. The exit condition of the loop was chosen when the main frequency’s signal to noise ratio went under 4. We considered that below this limit, all that is left is noise.

The process has been applied for each of the 14 blocks of continuous observations, leading to an average number of 43 different pulsations for each block. After identifying all these pulsations, we subtracted them from the photometric measurements. The resulting residuals were then rebinned from an initial time sampling of 120 seconds to a time sampling of 1800 seconds.

To validate the result, we checked that the three exocomet transits already identified in Ref.[20] are clearly visible in the data set cleaned from the Scuti variations by this procedure (Extended Data Fig. 2). The data are clean enough that new potential exocomet transits are suspected from the resulting light curve.

Remaining slow variations

After removal of the Scuti variations, the resulting light curve still shows slow variations. These variations can have various origins (residuals of the Scuti variations that were not properly eliminated, systematic correlated noise, instrumental effects, etc.). Whatever the origin of these slow variations, astrophysical or instrumental, they need to be cleaned before searching for exocomet transits. To do so, for each epoch of continuous observations we searched for a smooth function modelling these variations. First we rebinned the time series at a 1-day-long interval, in order to “protect” the shorter dimmings (such as comets), excluded the 1-day measurements where the deepest exocometary transits are already clearly identified (by using the detection procedure described in the main text), then we interpolated a function with a cubic spline to model the generic form of the time series. Finally, we normalized our time series with this model. The result is a flattened light curve, which can be directly used to search for shallow dip events due to exocomet transits. The procedure is illustrated in Extended Data Fig. 3 where a potential second shallower exocomet transit closely follows another deep exocomet transit event.

Comet’s nuclei radius

The conclusion on the size distribution given in the main text does not require to estimate the true size of the comet’s nuclei. It only relies on the assumption that the absorption depth is a good proxy for the dust production rate, which is supported by numerical simulations of exocomets light curves[18], and that the production rate is proportional to the area of the comet’s nucleus, which is consistent with Solar system comets’ models and observations.

Nonetheless, using these simulations and Solar system observations, we can make a step further and derive the typical sizes of the comets detected in the  Pic TESS light curve. From a newly calculated library of exocomet transit light curves similar to the one of Ref.[18], we derive a typical scaling law for the absorption depth AD, which iswhere is the dust production rate of the comet when it is at 1 au from the star, q is the orbital periastron distance and is the mass of the star. Because the transiting comet is seen from a distance that is extremely large compared to the star-comet distance, for a given cometary tail size, thickness and impact parameter, the percentage of the stellar disk covered by the tail and therefore the absorption depth do not depend on the comet distance to the star. The above estimate of the absorption depth is valid over a wide range of impact parameter and wide range of longitude of perisastron of the comet’s orbit, .

The periastron distance of the detected comets can be estimated using the transit time . We have , where is the length of the transit chord of the planetary trajectory on the stellar disk that is crossed by the planet at the velocity . The chord length depends on the transit impact parameter : it is twice the planet radius for a central transit with a zero impact parameter and it is zero for a tangential transit with an impact parameter equals to the stellar radius. Here we used the mean chord length obtained assuming a random uniform distribution of the impact parameter from zero to one stellar radius, that is . The transit velocity is given by . With a  Pic radius of (Ref.[41]) and a mass of , we find ) hours. The best fits values of correspond to distances ranging from 0.03 to 1.3 au, in good agreement for the distances expected for the comet evaporation. The mean value of the estimated periastron distances is about 0.18 au. Using this mean distance, we obtain the following relationship for the observed absorption depth and the dust production rate :

Finally, the relation between the evaporation rate and the comet’s nucleus size can be derived by scaling the observation of the Hale-Bopp comet. Using a radius of about 30 kilometers[23,42,43] and a dust production rate of at 1 au (Ref.[23]) for this well-observed dusty comet, we find . With a  Pic luminosity of , we find

All together, we conclude that the radius of the  Pic comets nuclei can be estimated using the photometric transit absorption depth with

Using this relationship, we derive a size of 1.5 km for the smallest detected comets (), and 6.7 km for the largest comet (). These sizes are remarkably similar to the sizes of comets in the Solar system.

Supplementary Information 1.