Exocomets size distribution in the β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} Pictoris planetary system

The star β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} 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 β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} 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 dN(AD)∝AD-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dN(AD) \propto AD^{-\alpha }$$\end{document}, where α=2.3±0.4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =2.3\pm 0.4$$\end{document}. This distribution of absorption depth leads to a differential comet size distribution proportional to R-γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{-\gamma }$$\end{document}, where γ=3.6±0.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =3.6 \pm 0.8$$\end{document}, showing a striking similarity to the size distribution of comets in the Solar system and the distribution of a collisionally relaxed population (γD=3.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _{{\text{D}}}= 3.5$$\end{document}).


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 Scientific Reports | (2022) 12:5855 | https://doi.org/10.1038/s41598-022-09021-2 www.nature.com/scientificreports/ 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 v sin i 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, t , 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 t 0 , 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 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 t 0 and using various plausible values for the t and β parameters characterizing the shapes of the exocomet transit light curves. We used β from 2 to 20 days −1 and t 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 t 0 . 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 χ 2 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 AD = K(1 − exp(−β�t)) . 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 δt as a function of the absorption depth can be fitted by a power law in the form dN(AD) = N 0 · (AD/10 −4 ) −α · (δt/100 days) · (dAD/10 −4 ) . Considering 29 events detected in 156 days of observations, that is all the 30 events except the deepest one (see below), we find that α = 2.3 ± 0.4 and N 0 = 33 +16 −11 , where the uncertainties have been evaluated using a Poisson distribution for the number of events in each bin of width dAD = 1.5 · 10 −4 .
The deepest transit event 20 of Julian Day JD = 2457000 + 1486 looks exceptional with an absorption depth of about 20 × 10 −4 . 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][11][12][13][14][15][16][17][18][19][20] www.nature.com/scientificreports/ 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 R 2 , we find that the differential number of exocomets as a function of the nucleus size is given by 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 γ = 3.6 ± 0.8 . 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 ( γ D = 3.5).
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  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. www.nature.com/scientificreports/ 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 γ = 3.6 ± 0.4 (stat.) ±0.7 (sys.), which is in remarkable agreement with our value for β Pic exocomets. For smaller long-period comets, a shallower distribution is found with γ = 0.35 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 γ = 3.7 ± 0.3 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 γ = 2.9 ± 0.2 for nuclei with radius larger than 1.25 kilometers, and γ = 1.2 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 γ = 3.5 power-law but truncated at small nucleus radii below 2.0 kilometers. In Ref. 29 a similar distribution is obtained with γ = 2.9 for radius between 2 and 5 kilometers, interpreting the distribution shallower than the canonical Dohnanyi's size distribution γ D = 3.5 (Ref. 24 ) as due to fragmentation of the JFC objects. In Ref. 30 the measured γ = 3.3 ± 0.2 for Jupiter family comets of radius between 2 and 10 kilometers is to be compared to the γ = 2.0 ± 0.1 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 γ = 3.55 ± 0.04 for the full sample with radius between 2.8 and 7 kilometers, γ = 3.2 ± 0.04 for near Earth objects (NEO) with radius down to 1.4 kilometers and γ = 3.45 ± 0.04 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 ( γ D = 3.5 ), 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 10 −3 , 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 Cleaning the light curves. The δ Scuti variations. β Pictoris is prone to δ Scuti type photometric variations [38][39][40] . These variations have a dominant frequency of 47.44 d −1 (corresponding to a period of about 0.5 hours), with an amplitude of up to 4 × 10 −3 (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 100 day −1 , 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 is where Ṁ 1 au is the dust production rate of the comet when it is at 1 au from the star, q is the orbital periastron distance and M * 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 t . We have �t = L chord /v transit , where L chord is the length of the transit chord of the planetary trajectory on the stellar disk that is crossed by the planet at the velocity v transit . 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 L chord = πR * /2 . The transit velocity is given by v transit = GM * /2q(cos ω + 1) ∼ GM * /q . With a β Pic radius of R * = 1.7R ⊙ (Ref. 41 ) and a mass of M * = 1.75M ⊙ , we find �t ≃ 13( q/1 au) hours. The best fits values of t 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 : www.nature.com/scientificreports/ 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 2 · 10 6 kg s −1 at 1 au (Ref. 23 ) for this well-observed dusty comet, we find Ṁ 1 au ≃ 2 · 10 6 kg s −1 (R/30 km) 2 (L * /L ⊙ ) . With a β Pic luminosity of 8.7L ⊙ , 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 ( AD ≃ 10 −4 ), and 6.7 km for the largest comet ( AD ≃ 20 · 10 −4 ). These sizes are remarkably similar to the sizes of comets in the Solar system.

Data availability
The observational data used in this work are publicly available in the Mikulski Archive for Space Telescope (MAST). The data in the tables and the final cleaned light curves are publicly available on GitHub at https:// github. com/ lecav eli/ BetaP ic_ TESS. www.nature.com/scientificreports/