Disk dissipation, giant planet formation and star-formation-rate fluctuations in the past three-million-year history of Gould's Belt

Although episodic star formation (SF) had been suggested for nearby SF regions, a panoramic view to the latest episodic SF history in the solar neighborhood is still missing. By uniformly constraining the slope $\alpha$ of infrared spectral energy distributions (SEDs) of young stellar objects (YSOs) in the 13 largest Gould's Belt (GB) protoclusters surveyed by Spitzer Space Telescope, we have constructed a cluster-averaged histogram of $\alpha$ representing YSO evolution lifetime as a function of the $\alpha$ value. Complementary to the traditional SED classification scheme (0, I, F, II, III) that is based on different $\alpha$ values, a staging scheme (A,B,C,D,E) of SED evolution is advised on the basis of the $\alpha$ statistical features that can be better matched to the physical stages of disk dissipation and giant planet formation. This has also allowed us to unravel the fluctuations of star formation rate (SFR) in the past three-million-year (3 Myr) history of these GB protoclusters. Diverse evolutionary patterns such as single peak, double peaks and on-going acceleration of SFR are revealed. The SFR fluctuations are between $20\%\sim60\%$ ($\sim40\%$ on average) and no dependence on the average SFR or the number of SFR episodes is found. However, spatially close protoclusters tend to share similar SFR fluctuation trends, indicating that the driving force of the fluctuations should be at size scales beyond the typical cluster sizes of several parsec.


INTRODUCTION
Most stars formed in clusters (Lada & Lada 2003;Porras et al. 2003), but it is still unclear how star clusters form in molecular clouds. Recently there has been increasing evidence that star formation in clusters is episodic. For example, multiple episodes of star formation or apparent acceleration of star formation have ucus (Alves de Oliveira et al. 2012), Perseus (Azimlu et al. 2015) and Taurus (Luhman et al. 2010), which indicates diverse star formation (SF) histories in these regions. Hsieh & Lai (2013) even briefly compared the α histograms of several nearby clouds. However, these α peaks are usually either not discussed in depth or simply attributed to SF episodes. The important fact that disk lifetime can be very different at different α values is usually neglected. Improved astrometric measurements provided by the GAIA satellite have revealed multiple episodes of SF also in older populations of some nearby SF regions such as Aquila/Serpens (Herczeg et al. 2019), Lupus (Galli et al. 2020), and the newly identified Orion C (Kounkel & et al. 2018). It is thus timely to exploit the power of low mass stars from large surveys (Megeath et al. 2022) to obtain a panoramic view to the episodic SF history in the earliest stage in the nearby protocluster forming regions.
The infrared SEDs of YSOs were classified using their slopes α by Lada (1987); Adams et al. (1987) and later improved by André et al. (1993); Greene et al. (1994) into Class 0, i, f, ii, and iii, which are roughly linked to the evolutionary sequence of YSOs. The α slope is a good tracer of the evolutionary time scales of circumstellar disks. Wilking et al. (1989) tried to assign Kelvin-Helmholtz timescales for the different SED classes, particularly the 2 Myr lifetime for Class ii sources. The life times of the SED classes were further refined by Evans et al. (2009) using the Spitzer c2d survey (PI: N. J. Evans) of nearby protostars.
However, our unpublished preliminary analysis of the star number counts of the different SED classes in the GB protoclusters has shown limitations in using the SED classification scheme to trace SF history. For example, the fractions of Class f YSOs seem to be correlated with that of Class 0+i ones among the GB protoclusters, which means the Class f is perhaps not a distinctive evolutionary stage of YSOs than Class 0+i. This point will be also confirmed by the correlation analysis in Appendix B. In this work, such limitations will be overcome by directly using the statistics of α, instead of the traditional SED classification, in the analysis of YSOs evolution.
In this paper, we will introduce the GB protoclusters from literature and the reprocessing of their IR SEDs in Sect. 2, analyze the α histograms, convert them into disk-evolutionary-age histograms, discuss the obtained disk-age histograms and average α histogram in Sects. 3.1 to 3.4, define the SED evolution stages and match them to disk dissipation and giant planet formation stages in Sect. 3.5, and finally summarize the main findings in Sect. 4.

Star samples and IR data
In order to understand the protoclusters in Gould's Belt as a whole, we collect published IR photometry between 2 and 24 µm of all GB protoclusters in the catalog papers, Rebull et al. (2010); Megeath et al. (2012) of Spitzer Taurus and Orion surveys and Dunham et al. (2015) of the Spitzer c2d (PI: N. J. Evans) and Spitzer GB (PI: L. E. Allen) surveys. These works have utilized color and magnitude criteria and, in some cases, spectroscopy to identify YSOs, excluded most background galaxies, outflow knots and part of field giant stars that also show IR excess, and matched the sources to 2MASS catalogue (Skrutskie et al. 2006). Several protoclusters in these papers are composed of more than one spatially well separated parts. We treat each part as a single protocluster and use their orientations (e.g., E for East, W for West, NE for Northeast) to differentiate them.
However, the contamination of giant stars to the Class iii YSOs is still an issue in these catalogs. We thus use the GAIA distances to identify and remove potential giant star contaminants as much as possible. For this purpose, all the above Spitzer YSOs are checked in the CDS X-Match Service 1 (Pineau 2020) to search for their GAIA counterparts in the Early Data Release 3 (EDR3) catalog (Lindegren et al. 2021) . A matching radius of 2 is used. We have found 4181 one-toone matches among the total of 6443 candidate member stars. For 299 multiple-matching cases, we select the closest GAIA source for the YSOs. Usually, mainly Class ii or iii YSOs are found to have GAIA counterparts. For the remaining 1963 candidate stars that have no GAIA counterparts, they are very likely embedded sources in the dense clouds so that we do not have accurate distance for them; we simply adopt them as member stars.
We adopt the following conservative criterion to identify nonmembers of each protocluster. First we observe the GAIA distance distribution of each protocluster by eyes and manually determine a reasonable distance range for member stars. In most cases, member stars distribute in a single continuous range of distances. However, in the few protoclusters near to Galactic plane (Aquila-Main, Serpens-Main and Serpens-NE), the candidate member stars gather around two distinct distance ranges. In these cases, we select the closer distance range for the protocluster, because the farther one is very likely dominated by background giant stars (contaminants). The adopted distance ranges are listed in Appendix Table 1. If a star's distance, including the corresponding distance uncertainties, lies entirely outside this distance range, then it is identified as a nonmember and excluded from our analysis. As a result, we have excluded 3 ∼ 354 contaminants from each of the 13 largest (with ≥ 50 final member stars) GB protoclusters (the detailed numbers of foreground and background contaminants are also give in Appendix Table 1). The other smaller GB protoclusters will not be considered further in this work. We have noticed that the adopted distance ranges are usually much larger than the transverse sizes (several pc) of the corresponding protoclusters, which we temporarily attribute to the large uncertainties of GAIA distances. We have examined the GAIA proper motions of the selected member stars and confirmed that the majority of them always concentrate well in a single small region.
The catalog paper of Dunham et al. (2015) also provided the extinction corrected SED slope α for the cluster member stars. However, as we will discuss in next section, the derived α histograms showed bewilderingly diverse peak positions, which is hard to interpret from a theoretical point of view. Thus, we decide to reconstruct the SEDs starting from the IR fluxes in the three catalog papers, using the extinction law of Weingartner & Draine (2001), and compute the SED slope α in a uniform manner for stars from all the three papers.

Extinction correction
The catalog papers of the GB protoclusters have estimated foreground extinction for each member star. In each SF region, the extinction was mainly determined for sources of SED Class ii and iii and then a uniform average extinction from these sources was adopted for embedded sources of Class 0, i and f. Therefore, the extinction in the literature mainly reflects the extinction of foreground dust. We adopt the extinction estimates A V or A K from the catalog papers, but use the extinction curve of Weingartner & Draine (2001, ; the extinction curve for R V = 5.5 for molecular clouds and Case B grain size distribution are obtained from a website 1 ) to convert A V or A K to the extinction in the other IR bands. Therefore, we need not to make the correction for scattering separately, as did by Dunham et al. (2015). We caution that the local extinction of embedded sources still remains after this step.

Saturation issues
1 The model extinction curve is from: https://www.astro.princeton.edu/ ∼ draine/dust/dustmix.html Most of the GB protoclusters have the Ks band and 24 µm photometry available for almost all their member stars, with the exception of Orion A, B-M and B-S due to a saturation issue in the Spitzer /MIPS maps. We have tried to recover the saturated 24 µm fluxes from their empirical correlations with the fluxes at shorter wavelengths among the unsaturated sources for 30%, 21% and 43% of member stars of the above three protoclusters respectively (see the details in Appendix A). A portion of member stars (10%, 3% and 19%, respectively) still have to be abandoned because photometry is absent for all wavelengths λ ≥ 5.8 µm. In addition, we drop a small number of member stars whose 24 µm flux is missing in several other protoclusters; the numbers of dropped stars are 9 (Ophiuchus), 2 (Serpens-Main), 8 (Taurus) and 5 (Perseus-W), which does not impact the statistics.

Fit the SED slope
The slope α of an IR SED is derived by fitting a straight line to it in the log (λF λ ) vs log λ plot. We choose to fit the SED from the 2MASS Ks band (2.2 µm) to Spitzer MIPS 24 µm band because the former traces the hot dust in the inner disk and the latter traces the cold dust in the outer disk (and envelop in earliest stages) so that the derived α will be a good tracer of the evolution of all circumstellar matter. For 528 stars that have the K-band flux missing (nearly 80% of them are Class 0, i, f, likely due to extinction), their SEDs start from the IRAC 3.6µm band. This does not impact the α values much. In addition, following this traditional definition of α offers us the opportunity to use the classical lifetime of 2 Myr of Class ii YSOs (Wilking et al. 1989;Evans et al. 2009) to do an absolute calibration to the timescale in the α-age conversion relation (see details in Sect. 3.3).

Two versions of SED-slope histograms
We initially used the extinction corrected α values from Dunham et al. (2015) to make the histograms. However, the peak positions of the histograms turned out to be significantly different among the protoclusters (the dashed histograms in Fig. 1). We stress that, because the histogram bins are very broad for smaller clusters, the shift of the α peak by one bin is already very large. This is hard to understand from a theoretical point of view because, if the evolution of SED is dominated by the common disk physics, as expected, the α-histograms of all protoclusters should peak around the same α value that corresponds to the slowest evolutionary stage. Only extremely sharp star formation rate (SFR) fluctuations have the potential to significantly shift the peak position, which means that the largest peak shift should appear in the narrowest α histogram. However, this is not the case in Fig. 1. For example, Serpens NE has the largest peak shift with respect to the others (among the dashed line histograms) but has a very broad α peak.
In order to resolve this problem, we start from the IR fluxes in the literature and redo the extinction correction and SED fitting in a uniform way for all protoclusters, as described in Sect. 2.2. We compare the derived α-histograms with that made from the literature α values in Fig. 1 for the nine protoclusters that have ≥ 50 member stars in Dunham et al. (2015). The effects of our extinction correction to these stars are the decreases of α values by 0.08 ∼ 0.37 on average (see the numbers after the cluster names in the figure). The histograms of the literature α values (dashed lines) generally peak at smaller α values than ours (full lines) and these peak shifts are roughly comparable to or even larger than the average effects of the extinction correction itself for many of these protoclusters, which reflects the over-correction of the dust-scattering extinction by Dunham et al. (2015).

All the SED-slope histograms
Our version of α-histograms of all the 13 GB protoclusters, including those in Orion and Taurus SF regions, are shown together in the left panel of Fig. 2. The most salient feature in the figure is that all the protoclusters show the major peak around α ∼ −1, indicating that these histograms mainly depend on the disk lifetime in different evolutionary stages traced by α. However, individual protoclusters also show secondary diversities in the histograms. For example, the major-peak positions can vary by ∆α ∼ 0.2 from cluster to cluster and several protoclusters show a sharp peak (e.g., Chamaeleon i, IC 5146-E, Orion B-S and Perseus-E), hinting on bursty SF events. A few protoclusters also possess older Class iii member stars with α < −2 (e.g., Serpens NE, Aquila-Main, Ophiuchus and Taurus) which means that their recent star formation might be initiated earlier than the other protoclusters.

Conversion between SED slope and disk age
The good agreement of the α histrograms among all the 13 GB protoclusters allows us to derive the disk evolution lifetime in the range of each bin of the histograms by averaging the histograms (after binned to the same α grid with bin width ∆α = 0.5) among the clusters to suppress the minor diversities. The so derived average α histogram, after a spline interpolation to a finer α grid  (Dunham et al. 2015, ;dashed lines) and from this work (full lines) for the nine protoclusters with ≥ 50 member stars in that paper. All histograms are normalized to an area of unity and vertically shifted for clarity. The number after each cluster name is the average amount of decrease of the α values of all member stars of the corresponding protocluster due to our extinction correction. The vertical gray full line marks the average α peak position of our version of histograms to facilitate comparison.
to facilitate later analysis, is shown as the black curves in the left half of Fig. 2 and the values of the average α histogram are listed in Appendix Table 2 to facilitate its reuse in future works. The difference of individual histograms with respect to the average can be interpreted as due to the temporal variation of SFR of each individual protocluster. The peaking of α ∼ −1 was found by other authors in individual protoclusters such as Auri-gaCMC (Broekhoven-Fiene et al. 2014), Cepheus (Kirk et al. 2009), Chamaeleon i (Luhman et al. 2008), Ophiuchus (Alves de Oliveira et al. 2012), Perseus (Azimlu et al. 2015), and Taurus (Luhman et al. 2010), but lit- Figure 2. Histograms of spectral energy distribution (SED) slope α (Left) and disk age (Right) of YSOs of the 13 largest protoclusters in Gould's Belt as surveyed by Spitzer Space Telescope. The number of member stars is shown behind each object name. The black curve in the left panel is a smoothed version of the average histogram of the 13 protoclusters. The averaging was only performed to the right side of the vertical black dashed line, i.e. within −2 ≤ α ≤ 5, and the common bin width for averaging is ∆α = 0.5. All histograms are normalized to an area of unity in the considered value range and vertically shifted for clarity. The vertical black full line marks the theoretical limit of α = −3 for main sequence stars much hotter than 3000 K. The disk-evolution-time histograms in the right panel are converted from the α histograms according to Sect. 3.3. Negative time means in the past and 0 Myr is now. The error bars are the dispersion among the 13 protoclusters for the black curve in the left panel and the square root of star numbers in each time bin in the right panel. Both the short vertical red lines on the α-axis in the left panel and the long vertical red dashed lines in the right panel mark the border lines of traditional classification of protostar SED types: 0+i, f, ii and iii. The bin width of each histogram is optimized with the 'Freedman-Diaconis' rule (Freedman & Diaconis 1981).
tle attention was paid to test its link to disk evolution lifetime.
However, there are several weaknesses in the above treatment. First, if all the protoclusters somehow share the same temporal variation trend in α during their formation, this common trend will enter the averaged his-togram and alter the derived disk evolution lifetimes. Second, the number of protoclusters is still small so that uncertainties are still large in the average α histogram. Thirdly, there exist SED diversities either due to different disk evolution paths (Strom et al. 1989;Lada et al. 2006;Cieza et al. 2007) or due to random disk inclina-tion angles Crapsi et al. 2008). Finally, the number of the youngest member stars (with α > 1) could be underestimated due to large extinction to their locations in the dense clouds.
Despite these weaknesses, it is still instructive to integrate the distribution of disk-evolution lifetimes (represented by the α-bin heights) to transform the α histograms into that of the evolution age of an "average disk". To integrate the average histogram from α = 5 to any α value between −2 and 5, we first fit the histogram with a second order spline function using the splrep and splint functions of the scipy.interpolate B-spline package (Note: A cubic spline will produce an unrealistic oscillation). The spline fit is denoted as P (α) and normalized into a probability density function as 5 −2 P (α)dα = 1. Then, the disk life time in any unit α range of dα is dt = aP (α)dα, where a is a constant scaling factor to be calibrated. The conversion relation from α to disk evolution time T (or negative disk age) can be numerically derived as A graphic presentation of Eq. 1 is given in Fig. 3 and the data of this curve are given in Appendix Table 2 for reuse in future works. The error bars are numerically evaluated by a Monte Carlo procedure in Appendix B. We can see that, for an average disk, the α value first decreases quickly with time in the early stages when α ≥ 2.5. Then, the decrease of α gradually slows down, reaches a slowest point around α = −1 and speeds up again towards α = −2. These trends have fundamental indications to the dissipation rate of circumstellar dust at a population level, which we will discuss in detail in Sect. 3.5. We stress that this conversion relation is less meaningful for individual stars than for star populations because of SED diversities, as we will discuss later.
Applying Eq. 1 to every member star of all the GB protoclusters, we derive the histograms of disk evolutionary ages in the right half of Fig. 2. The error bars of the disk-age histograms come from the square root of star number of each bin. We do not transfer the error bars in Fig. 3 into the disk-age histograms, because these error bars reflect both the star-number counting error of each α histogram and the error due to the fluctuations of SFR, while the purpose of the right half of Fig. 2 is merely to present the same SFR fluctuations through Figure 3. The conversion relation from SED slope α to disk evolutionary age for YSOs, as expressed by Eq. 1. The error bars are numerically obtained as in Appendix B. The data behind the plot is given in Appendix Table 2. directly comparing the disk-age histograms. However, if these disk-age histograms would be compared with similar histograms obtained with other methods in future, the error bars in Fig. 3 should be transferred to the disk-age histograms through another simple Monte Carlo approach to replace the current number counting error bars.
We comment that the derived disk-age histograms look very different from the original α histograms in the left half of Fig. 2. This is because, on the one hand, the corresponding ranges of disk-age are distinct in different α bins and, on the other hand, the disk-age histograms only reflect the differences between the individual α histograms and the averaged one.
We also remark that our results indicate apparent disk-age uncertainties better than 0.5 Myr in the youngest stage of < 3 Myr, as can be seen through the error bars in Fig. 3. However, we have implicitly assumed that all YSOs have the same disk-evolution track traced by α. If other factors such as the different stellar mass, disk mass, accretion history, metallicity, initial angular momentum, mass reservoir, etc. would be taken into account, the absolute disk-age uncertainties will be larger. The situation is similar to stellar ages constrained through the Hertzprung-Russell Diagram (HRD) for similarly young stars where uncertainties in the theoretical models make the stellar age less accurate than apparently determined on the HRD (Soderblom et al. 2014). As we will discuss in a later subsection, some of these neglected factors are responsible for the observed SED diversities.
We caution again that the sample of the youngest and reddest member stars are usually embedded in dense cloud cores and are thus incomplete due to high local ex-tinction. Thus, the real evolution of the SED slope during the youngest stages could be slower than we found here.

Disk-age histograms
The derived disk-age histograms in the right half of Fig. 2 clearly demonstrate, for the first time, the acceleration and/or deceleration of SFR in the latest 3 Myr history of each protocluster. Some protoclusters (Serpens NE, Serpens-Main, Aquila-Main, Ophiuchus and Chamaeleon i) are older, experienced their peak SF epoch more than 2 Myr ago. Some younger protoclusters (Perseus-E, IC 5146-E, Orion A, B-M, B-S and the major peak of Taurus) reached the maximum SFRs between 1 ∼ 2 Myr ago. The two youngest protoclusters, Perseus-W and AurigaCMC, have their SFRs accelerating up to now. Some protoclusters seem to show multiple peaks in the disk-age histograms. A dip test of Hartigan & Hartigan (1985) 2 is performed to the diskages of each protocluster. If the p-value is < 0.05, the tested sample is multi-modal, otherwise if the p-value is < 0.1, the tested sample possibly has multiple modes. Aquila-Main (p = 0.0), Ophiuchus (p = 0.004) and Taurus (p = 0.008) are found to be multi-modal and Serpens-Main (p = 0.09) possibly multi-modal, indicating multiple SF episodes even within the short period of 3 Myr. Combining these results with episodic SF events discovered on longer timescales in older Pre-Main Sequence (PMS) populations selected from other surveys than Spitzer (e.g., Beccari et al. 2018;Herczeg et al. 2019) highlights a picture of multi-timescale episodic star formation.
An interesting feature in Fig. 2 is that spatially close protoclusters tend to share similar SFR flucturation trends. For example, Serpens NE, Aquila-Main, Serpens-Main and Ophiuchus protoclusters are close to each other on the sky plane and all of them reached the SFR peak early (about 3 Myr ago). The Orion A, B-S and B-M protoclusters are also close neighbors and their SFRs peak around the similar time of ∼1.5 Myr ago. The AurigaCMC and Perseus-W protoclusters are also not far from each other and both of them are experiencing the acceleration of SF even up to today. The two protoclusters with the most prominent double peaks in the disk-age histograms, Aquila-Main and Ophiuchus, are also in adjacent regions. These features indicates that the fluctuations of SFR should be driven by agents at large size scales covering several adjacent protoclusters.
However, we also note that spatially close protoclusters do not necessarily always have similar SFR fluctuation trends. For example, Perseus E and W show the SFR peak at quite different times. This is also consistent with the diversity of SFR trends among sub-clusters within the same protoclusters, which we will discuss in our next paper. It hints on the complexity and diversity of cluster formation.
In order to quantify the fluctuations of the SFRs ( ), we compute the standard deviation of each disk-age histogram and show their ratios to the average height of the corresponding histograms in Fig. 4. Because the diskage histograms are normalized to total areas of unity, the height of each bin is proportional to the level of SFR in each disk-age bin. Thus, the standard-deviation-toaverage ratio of a disk-age histogram is also the ratio of the standard deviation (σ ) of its SFRs to its average SFR (¯ ) or, in another word, its relative fluctuation of SFR (σ /¯ ). One can convert the SFRs into a unit of M yr −1 by multiplying them with an average stellar mass of 0.5 M . It can be seen that the SFRs fluctuate strongly with amplitudes between 20% to 60% (∼ 40% on average) and the fluctuations are not correlated with the average SFRs, nor with the number of peaks in the disk-age histograms. On the other hand, despite the large SFR fluctuations, all protoclusters maintain a non-negligible level of SF most of the time so that the SF should be a continuous activity within timescales of 3 Myr, instead of neat discrete bursts.

A uniform statistical view of disk dissipation and giant planet formation
In this section, we will first show that the shape of the average α histogram has not been impacted much by the sample selection biases. Then, we will introduce a new SED staging scheme that has the potential to neatly match the statistical features in the average α histogram to individual physical stages of disk dissipation. Finally, we will propose a scenario to naturally integrate the giant planet formation into the new scheme to play a key role in the disk dissipation.

Sample biases
For the first time, we have been able to split the effects of disk evolution timescales and fluctuations of SFR in the α histograms of protoclusters at a population level. The common trends in the α histograms in Fig. 2 or the α-age conversion curve in Fig. 3 shall cast new light on the dissipation processes of circumstellar dust around YSOs. However, before discussing that, we need to clarify how these trends have been distorted by selection effects. Beside the possible incompleteness of the youngest member stars (with α > 0) due to high (uncorrected) opacity in the dense molecular clouds, as mentioned in Sect. 2.2, we will discuss another two biases below.
The first possible bias is that some YSOs might have been removed together with the background galaxies. Compared to typical YSOs, galaxies usually show redder colors at a given 24 µm flux level (e.g., Harvey et al. 2007). Thus, the confusion with galaxies might have enhanced the incompleteness of the reddest YSOs of Class 0 and i to some extent. However, because the YSOs and galaxies are separated well on the multiple color-magnitude plots, this bias should be small.
The second bias is related to the removal of main sequence (MS) stars when selecting YSO candidates from the original Spitzer maps. For example, when Harvey et al. (2007) selected YSOs in the Serpens regions, about 94% of the Spitzer point sources with sufficient fluxes to make informative SEDs had been removed as MS stars. According to their approach (described by Evans et al. 2003), all good fits with a theoretical photosphere model and the extinction law of Weingartner & Draine (2001) are considered as MS stars. This means the removed MS stars should always have their de-reddened IR SEDs very similar to a blackbody function whose IR SED slope is theoretically ∼ −3 and this bias mainly impacts this α range. However, the fact that the Taurus protocluster in our sample has shown a secondary α peak just around α = −3 in the left half of Fig. 2 indicates that this bias is not necessarily very serious even in this α range. Therefore, because our discussions will be focused in the α ≥ −2 range, we need not worry about this bias too much.

Stages of SED evolution and disk dissipation
It is thus reasonable to assume that the average α histogram is roughly free of observation biases and a physical interpretation is warranted. To facilitate discussions, we re-plot the smoothed average α histogram and the αage conversion curve together in Fig 5. Because the α value is mainly determined by the K-band and/or 3.6 µm emission from the hot inner disk (or from photosphere in the case of transition disk) and the 24 µm emission from the outer cooler disk beyond 5 ∼ 10 AU radius or even the infalling envelope in early stages, these curves mainly reflect the evolution of the entire protoplanetary disk and envelope from the largest positive α (∼ 5 in our sample) to the photospheric value α ≈ −3.
We can identify five α sub-regions that show distinct evolutionary characteristics in Fig. 5 (we call it SED staging scheme): Stage A (2.5 < α) shows a very rapid decrease of α with time at the earliest stage so that very few sources are found in this stage (very short timescale of ∼ 0.01 Myr) which should be true even after taking into account the sample incompleteness; Stage B (0 < α ≤ 2.5) witnesses a gradual increase of star number towards the transition region from the Class i to f (i.e., a gradual increase of disk dissipation timescale with time or a gradual slowdown of the decrease of α) in a period of ∼ 0.5 Myr; Stage C (−1 < α ≤ 0) covers the ascending half of the α-histogram peak which signals a rapid increase of star number toward the middle of the Class ii (i.e., an rapid increase of disk evolution timescale over time or a significant slowdown of the disk dissipation) in a period of ∼ 1 Myr; Stage D (−2 < α ≤ −1) strides across the Class ii and iii, oppositely containing the descending half of the α-histogram peak which indicates a sharp decline of star number (i.e., a sharp decrease of disk evolution timescale or a fast speedup of the dissipation of the remaining disk until its IR excess almost disappears) within another period of ∼ 1 Myr; Stage E (α ≤ −2) brackets the diskless PMS stars and fainter debris disks (the timescale is difficult to estimate due to the lack of observed star samples and the potential confusion with MS stars, as discussed above. The five SED Stages above very likely correspond to distinctive dissipation processes of the circumstellar dust at the level of YSO population. According to the current theory, a newly formed single star first appears as a Class 0 protostar that has its youngest disk still embedded in a gas envelope (André et al. 1993). The envelope will be lost during Class i stage (Lada 1987). The disk starts its long journey of dissipation mainly through the so-called UV-switch mechanism (Clarke et al. 2001) in which the disk first loses mass through viscous accretion onto the star and then FUV, EUV and X-ray evaporation of the outer flaring disk will take over to switch off the mass supply to the inner accretion disk and leave behind an inner hole (a transitional disk; e.g., Alexander et al. 2006;Gorti & Hollenbach 2009). At the same time, dust grain settling towards the mid-plane of the disk (e.g., Dullemond & Dominik 2004) and its growth to big grains and planetesimals (e.g., D'Alessio et al. 2001) should proceed in parallel. The evolution of disks may have started quite early within the first 1 Myr (Furlan et al. 2009). Thus, the earlier part of the disk dissipation should have started during Class i while its later part roughly corresponds to Class ii (classical T Tauri stars; CTTS). Finally, the photoevaporation, together with other processes such as stellar accretion, grain growth and planet formation, will quickly dissipate the remaining primordial disk during Class iii stage and a CTTS rapidly evolves through a weak-line TTS stage (WTTS; see e.g., Alexander et al. 2006;Williams & Cieza 2011). The observed dissipation timescale of the entire primordial disk is a little bit controversial in literature, ranging from 4 ∼ 8 Myr (Haisch et al. 2001;Ribas et al. 2014) down to ∼ 1 Myr (Cieza et al. 2007) and mostly around 2 Myr (Spezzi et al. 2008;Mamajek 2009;Richert & et al. 2018). Further on, some planetary systems may develop a debris disk through collisions of planetesimals on longer timescales ≥ 10 Myr (Wyatt 2008).
However, it has been suspected that not all disks follow the unique path of evolution delineated above. This was mainly inferred from the diversities of IR SEDs and millimeter maps (e.g., Furlan et al. 2006;Cieza et al. 2007;Owen 2016). The diversity of SEDs is also found in our sample of YSOs. We present the median SEDs of several subgroups of YSOs near the border lines of the SED stages and inside Stages B, C and D in Fig. 6. To create the median SEDs, we first fit a second order polynomial to each SED, then sort them according to the concavity (the second derivative) of the polynomial, and finally sequentially divide them into 3 or 4 roughly equal subgroups to compute median SEDs. Because YSOs around the border of Stages A and B is few, we include all Stage A YSOs into this group. The SEDs of this group have no detection in 2MASS Ks band so that they have to be normalized to the 3.6µm band. The α = 1 ± 0.1 group also has some Ks-band non-dections and only the 36 member stars with Ks-band flux are considered. In addition, some member stars in the Orion region have no 24µm flux due to the MIPS saturation issue discussed in Sect. 2.3. They are omitted in the computation of the median SEDs, which mainly impacts Stage C (α = −0.5 ± 0.1 and −1 ± 0.1). We will briefly discuss the major trends below, while more details will be analyzed in a dedicated paper.
As can be seen among the median SEDs, the diversity has already been large from the very beginning in Stage A and begins to diminish from the end of Stage B (α = 0) and on. A low level of SED diversity is maintained throughout Stage D. The well known transitional disks discovered since the work of Strom et al. (1989); Skrutskie et al. (1990) should have the most concave SEDs in our sample, because they usually show little excessive emission at NIR and MIR (λ < 20 µm) but maintain significant excess beyond 20 µm as optically thick full disks do (Espaillat et al. 2014). They mainly appear in our Stages C and D and are all wrapped in the most concave median SED in the plots (the bottom red curves in the panels with α = −0.5 ± 0.1, −1.0 ± 0.1, and −1.5 ± 0.1). Thus, no more than 25% of the sample YSOs in any α sub-range of the two stages belong to transitional disks.
The SED diversity during Stage A can be attributed mainly to geometrical features of a standard 'enve-lope+disk' model. According to the SED evolutionary models of Whitney et al. (2003), the flaring disks in Figure 6. Representative median SEDs normalized to 2MASS Ks band (or IRAC 3.6 µm band when the Ks band is nondetection) around each border line of the SED stages in Fig. 5 and inside Stages B, C and D (see the selected α ranges and the number N of YSOs in each group in the figure legends). The method to define the median SEDs can be seen in Sect. 3.5.2. The sequence of panels are in the same order of α values as in Fig. 5. Blue, yellow, green and red colors are used sequentially for the median SEDs from the most convex to the most concave in each panel. The black dashed curve is a 6000 K black body.
the earliest evolutionary stage should be embedded in an opaque infalling envelope with bipolar cavities. The opening angle of the cavities and the viewing angle to the cavity axis have strong impact to the observed SED morphologies. The model SEDs seen in different viewing angles in their late Class 0 and Class i cases (their Fig. 3) are sufficient to explain the median SEDs we found.
The SED diversity during Stage B needs more factors to account for. The median SEDs show convex, flat and concave shapes. The latter two shapes can be similarly explained with the model SEDs viewed at different inclination angles in the Class i and late Class i cases of Whitney et al. (2003). However, the convex SEDs need either additional emission in the 3 ∼ 8 µm range or deficit of emission in Ks and/or 24µm bands to explain. This cannot be interpreted solely by the puff up of disk inner rim, because it mainly enhances emission around 2 ∼ 3 µm for massive disks but not for T Tauri stars (Dullemond et al. 2001). We speculate the following candidate mechanisms for further work to test: 1) these YSOs with convex SEDs are likely heavily embedded sources so that the uncorrected local extinction preferentially attenuates the light at shorter wavelengths; 2) they could be binaries that have truncated the outer disk and reduced the 24 µm flux (Williams & Cieza 2011). However, only the local extinction scenario has the potential to also explain the diminishing trend of SED diversity towards later stages.
Because the SED diversities in the early stages can be satisfactorily interpreted as above by different cavity opening angles and viewing angles, the issue of multiple disk evolutionary paths mainly exists in the later evolutionary stages. Wood et al. (2002); Lada et al. (2006) proposed a 'homologously depleted disk' or 'anemic disk' model for a large population of observed Class ii disk sources that always maintain a power law shape of SED during the evolution. Our Fig. 5 also shows that the largest portion of the Class ii and iii YSOs in GB do have power-law SED shapes. Such disks cannot be in-terpreted alone by grain settling and growth and stellar accretion, because the fast dust settling and coagulation must be balanced by replenishment through turbulence in the disk (Dullemond & Dominik 2005) and the accretion may be too slow in later stages to explain the observed fast disk dispersal. However, combining them with the self-shadowing effects of an inner wall (it reduces the far IR disk emission; Dullemond & Dominik 2004) and photoevaporation (Alexander et al. 2006) has the potential to explain the entire observed evolution sequence of anemic disks. On the other hand, a small number of millimeter-strong transition disks (Owen 2016) represent another distinct evolution pathway of some Class ii and iii objects which may need the formation of giant planets to clear out the inner disk within 2 ∼ 3 Myr (e.g., Rice et al. 2003;Varnière et al. 2006;Dodson-Robinson & Salyk 2011;Andrews et al. 2011;Gorti et al. 2011;Dong et al. 2015;Sallum et al. 2015). The existence of multiple pathways may pose a question on the short lifetime of transitional disks that was constrained by counting stars also in the 'anemic disk' channel. The stars' choice between the two evolutionary paths might depend on the details of planet formation, which we will discuss in the next subsection.
Referring to the above picture, our observed SED evolutionary stages can be matched to disk dissipation processes as follows. The very few sources observed in Stage A are likely stars with infalling envelope actively interacting with bipolar jets. This is a brief stage of ∼ 10 4 yrs for the newly born stars to emerge at NIR. The viewing angle of the bipolar cavities strongly impact the detectability of the stars and their SED diversity.
Detectable disk sources begins to build up during Stage B. In this ∼ 0.5 Myr process, it is natural to expect the dissipation of the remaining envelope by widening the bipolar cavities, grain settling toward the mid-plane and grain growth. The flattening of the disk results in lower surface temperature, which explains the observed quick decrease of α. The disk flattening slows down when the disk becomes flatter, giving rise to the observed gentle rise up of the α histogram with decreasing α. The SED diversity due to the different bipolar cavity opening angles and viewing angles also begins to diminish with disk flattening. We suggest that the creation of an inner-disk wall might have started at a certain moment because the strong stellar accretion should be powerful enough to generate energetic radiation to blow up the wall simultaneously.
After then, the SED evolution further slows down significantly in Stage C. This is because a stable flaring disk structure is established during this ∼ 1 Myr stage. The disk should have an inner wall that shadows and protects a large part of the flaring disk. The disk evolution is mainly driven by the slow stellar accretion through the wall. Dust grain settling and growth is balanced by turbulent grain replenishment. Planetesimals and planets begin to form in the mid-plane. This stage ends at a yet poorly understood critical point around the α = −1 peak in which the slowing trend of disk dissipation reverses to an accelerating trend (we will discuss the mechanism in next subsection). The sub-millimeter mapping of YSOs in the Taurus-Auriga SF region by Andrews & Williams (2005) has uncovered a sharp decrease of disk mass during the Class ii stage (their Fig. 15), which also supports the α peak to be the critical starting point of the final runaway disk gas dissipation.
Stage D should be the final stage of the primordial disk. The further dissipation of the entire disk speeds up, possibly due to the (partial) clearing of the inner disk. In the average α histogram, the profile shape of this stage is roughly left-right symmetrical to that of Stage C, indicating a nearly reciprocal timing of the two stages. The division of Stages C and D is thus very different from the so-called two-time-scale division between the long full disk stage (roughly Class ii) and the brief transitional disk stage (roughly Class iii; e.g., Alexander et al. 2014;Ercolano & Pascucci 2017). Therefore, a different physical interpretation is needed here. We propose that, due to some unknown processes at the α = −1 peak, the inner wall begins to collapse and the disk is quickly dispersed by the combination of stellar accretion and photoevaporation in the ∼ 1 Myr period of this stage. The other processes such as grain settling and growth only play minor roles due to the diminishing gas and dust density.
Finally, the YSOs become diskless stars or debris-disks in Stage E. We comment that some debris disks could return from Stage E to Stage D if their 24µm dust emission would be strong enough. However, identifying such objects in our sample is beyond the scope of this work.

Role of giant planet formation in disk dispersal
The paramount α = −1 histogram peak (in the middle of the traditional Class ii) is not only the slowest disk dissipation moment, but the critical point when the disk dispersal reverses from slowing down to speeding up. No mechanism is known for this reversal. Processes such as stellar accretion, dust settling and grain growth are not expected to show such a critical state. Photoevaporation does possess a critical point when the evaporation mass loss rate overtakes the accretion rate and an inner hole will be cleared out and the outer disk will be shed quickly. However, this critical point is expected only after several Myr of stellar accretion (e.g., Clarke et al. 2001;Alexander et al. 2006;Gorti & Hollenbach 2009), which is too long compared to the age of the α = −1 peak (1.76 Myr). In addition, this is the mechanism for the two-time-scale problem of disk dispersal, as mentioned in Sect. 3.5.2, which does not agree to the observed symmetry in the average SED histogram between Stages C and D.
The only hopeful driving process is the widely studied giant planet formation. The giant planets in the Solar system, Jupiter and Saturn, have a huge gas sphere so that they must accrete enough gas before the dissipation of the gas disk. Therefore, in our average α histogram, the giant planets of low mass stars very likely mature around the α = −1 peak. Ercolano & Pascucci (2017) reached a similar conclusion in their review of disk dispersal and planet formation and even suggested to call Class ii sources as planet forming disks, instead of protoplanetary disks. The capability of giant planets in clearing the inner disk (so as to speed up the dispersal of the outer disk) is also supported by observations. For example, Andrews et al. (2011) found that only dynamical clearing by giant planets (or brown dwarfs) can explain the cavities in 12 transition disks they mapped at sub-millimeter; Sallum et al. (2015) caught accreting planets in action in the multiple ring system of the transition disk source LkCa 15. Although the formation of giant planets alone cannot explain the disappearance of disks because of their too low masses (Pascucci & Tachibana 2010), it can play a key role in combination with photoevaporation (Rosotti et al. 2013).
In order to better understand the capability of giant planets in speeding up disk dispersal, let us recall the typical three step scenario of giant planet formation (Pollack et al. 1996). In the first step, a runaway accretion of cm-sized pebbles builds up an initial solid core above ∼ 10 M ⊕ to enable gas accretion (Lambrechts & Johansen 2012). The second step is the most time consuming: the gas accretion will quickly clear up the feeding zone of the protoplanet to greatly slow down fur-ther accretion of solids and the stable gas sphere cannot grow further by adding more gas. However, the accreting solid core will dynamically migrate radially to maintain a non-negligible accretion rate; it will achieve a critical cross-over mass threshold (i.e., the accreted planetesimal mass equals the accreted gas mass) within about 1 ∼ 3 Myr (Alibert et al. 2005;Cieza et al. 2012;Bitsch et al. 2015). Finally, a runaway gas accretion is triggered (because the accreted gas sphere starts to gravitationally collapse; Mizuno et al. 1978) to complete the giant planet formation rapidly.
The runaway gas accretion in the last step has the capability to open a gap in both the dust and gas disks (Crida et al. 2006;Fung et al. 2014) and to play a key role in reversing the slowing trend of the disk dissipation. Although the planet opened gaps may not totally block out the supply of gas and small dust grains to the inner accretion disk (e.g., Rice et al. 2006;Crida & Morbidelli 2007;Zhu et al. 2012), it can significantly reduce the gas density of the inner disk to speed up the dispersal of the outer disk by photoevaporation (Dodson-Robinson & Salyk 2011;Andrews et al. 2011;Rosotti et al. 2013). This speed-up commences around α = −1 according to our average α histogram.
However, the above giant planet formation scenario is only obviously applicable to transitional disks, but not to anemic disks that do not show the signature of inner disk hole in their SED morphologies. There is no readily available model to explain why these disks also speed up their dispersal around α = −1. We speculate that the giant planet formation might still work in this case, but the planet masses could be lower than in the transitional disks so that the shallower gaps and inner holes are not discernible in SEDs. However, the shallow gaps and holes still reduce the mass supply to the inner disk to speed up the photoevaporation. Of course, other processes such as grain settling and growth could be contributing too.
There have been some evidences supporting the planet scenario for anemic disks. On the one hand, Najita et al. (2007);Sicilia-Aguilar et al. (2008); Muzerolle et al. (2010) found that anemic disks tend to appear around lower mass stars with lower disk masses while transitional disks around more massive stars with higher disk masses, indicating that the former tend to form less massive giant planets than the latter. This is also consistent with the very low occurrence rates of massive giant planets in direct imaging observations (e.g., Bowler 2016). On the other hand, imaging studies have been providing acuter clues. The ALMA large program DSHARP (Andrews et al. 2018) have revealed that 16 out of the 20 nearby Class ii anemic disk targets show gaps and/or spirals related to planet dynamics. The ALMA archival study of Francis & van der Marel (2020) have also shown that many disks with gaps and inner holes at millimeter wavelengths have their SED shapes actually closer to anemic disks than to typical transitional disks (their Fig. 16). With more ALMA continuum or CO line maps, Long et al. (2018);Wölfer et al. (2022) further proposed that many of the disk sources (including transitional disks) have their disk gaps possibly opened by low mass planets. All of them hint on the ubiquitous role of planet formation also in anemic disks defined from SEDs.
If the ubiquitous role of planet formation in anemic disks would be further confirmed, there emerges a simple uniform picture in which only the rare most massive giant planets (or even stellar companions) can clear out the massive inner dust disks to leave behind the small number of millimeter bright transitional disks reviewed by Owen (2016), while the more common less massive planets grow in and help disperse the larger population of anemic disks.
Finally, we comment that the current theory of giant planet formation still has the room for improvements. For example, many millimeter-bright transitional disks with a large inner disk hole are found to show sizable gas accretion (Manara et al. 2014;Owen 2016), whilst giant planet formation and photoevaporation models constantly predict too many non-accreting transition disks (Rosotti et al. 2015;Ercolano et al. 2021). Further studies are desired to solve such conundrums and to explain the α = −1 peak identified here.

SUMMARY
By reprocessing in a uniform manner the IR SEDs of all member stars of Gould's Belt (GB) protoclusters surveyed by Spitzer Space Telescope, we have successfully identified the common features in the histograms of SED slopes (α) that represent the distribution of disk evolution life time over the α value. This has allowed us to convert α into disk evolution age to obtain the first panoramic view to the large fluctuations of SFRs in the GB protoclusters in the past 3 Myr history. We have found the following properties: (1) The star formation is roughly continuous in GB and the amplitudes of SFR fluctuation are 20% ∼ 60%, or ∼ 40% on average; (2) Spatially close protoclusters tend to share similar SFR fluctuation trends, indicating common driving forces for the fluctuations at large size scales; (3) The amplitudes of fluctuation are likely independent of the average SFR and the number of SF episodes.
The average α histogram is largely free of sample selection bias. It shows several statistic features that allow us to define five SED evolutionary stages: A, B, C, D and E, with their dividing points at α = 2.5, 0, -1 and -2 or at negative SED evolution ages of -0.01, -0.55, -1.76 and -2.84 Myr. It can be considered as a further development of the traditional SED classification scheme of 0,i,f,ii and iii that was based on different α values and on the theoretical meaning of the different SED classes.
The SEDs of individual YSOs show diversities. The diversities in the earlier Stages A, B and C can be interpreted by different opening angles and viewing angles of bipolar cavities, while the diversities during Stage D indicates two disk evolution pathways: anemic and transitional disks. Despite the diversities, the SED staging scheme provides us a convenient time frame of disk evolution to coordinate further studies. The SED staging scheme can be better matched to the dissipation stages of circumstellar matter as follows: (1) bipolar cavities begin to develop in the infalling envelope in the first 10 4 yr; (2) the infalling envelope is dissipated by stellar jets and winds and grain settling and growth, and possibly an inner wall is puffed up within ∼ 0.5 Myr; (3) a stable flaring disk structure protected by the shadowing of an inner wall is developed and stellar accretion in the inner disk, photoevaporation of the inner wall, grain settling, growth and planet formation in the disk mid-plane jointly drive its slow evolution in ∼ 1 Myr; (4) giant planets form and open gaps in the disk to starve the inner accretion disk, which triggers the acceleration of the final disk dispersal by photoevaporation in another ∼ 1 Myr; (5) the final stage of diskless star or debris disk is not well constrained in this work. This work has also opened up the opportunity to have a better look into the evolution of the cluster spatial structures with time. We will report it in the next paper.
We thank the science referee for helping us improving several important defects and pushing us to take into account the diversities of SED and disk evolution paths, which has significantly improved our data analysis and discussions.
We thank Prof. Diego Mardones for many useful discussions in the early stage of the work. ML jointly initiated the project, performed all data collection, processing, analysis, code writing and participated in the manuscript revision. JH jointly initiated the project, guided every step of the entire work, wrote and revised the manuscript. The other co-authors participated many discussions and significantly improved many aspects of the data analysis and science discussions.
JH and TL acknowledge the support from the National Natural Science Foundation of China ( The Spitzer /MIPS maps of the Orion regions suffer from saturation issues in a few brightest areas so that no reliable 24 µm (in some cases also no IRAC 5.8 and 8.0µm) photometry is available for a significant fraction of member stars in those protoclusters (Orion A, B-S and B-M; Megeath et al. 2012). To partially remedy the bias resulted from the missing of these member stars, we need an approach to reconstruct their α(2.2 ∼ 24 µm) values from available data at shorter IR wavelengths. In each of the three protoclusters, we first compare the histograms of α values defined in 2.2 ∼ 5.8 and 2.2 ∼ 8.0 µm wavelength ranges between the member stars inside and outside of the saturated areas and find that they are statistically similar. It is thus reasonable to assume that the relationship of the α values defined in the three wavelength ranges are statistically similar for all member stars, no matter whether they suffer from the MIPS saturation issue. Then, for those member stars outside of the saturated areas, we compare their α values defined in different wavelength ranges and find the following good linear correlations in all cases.

B. ERROR BARS OF THE SED-SLOPE-TO-DISK-AGE CONVERSION CURVE
The error bars of the SED slope (α)-disk age conversion curve in Fig. 3 cannot be analytically obtained from the error bars of the average α-histogram, because it involves integration of the latter. We obtain them by numerically re-sampling the average α-histogram to imitate its error bars.
First, we compute the correlation matrix of the 14 α-bins of the average α-histogram from the original α-histograms of the 13 protoclusters. This step is necessary because the different α-bins are not necessarily independent when being averaged over the 13 protoclusters. The resulting correlation coefficients are shown in Appendix Table 3. The only valuable comment we can make is that the bins 5, 6, 7, 8, 9, 10 and 12 (α = 0.25 ∼ 3.75) are significantly positively correlated. This could be explained by the short disk life times in these bins so that these bins can be under the control of a single SF episode. This fact confirms our preliminary finding of the positive correlation of star number counting between Class f and 0+i protostars (not published).
Then, we generate 10 4 random realizations of the average α-histogram using this observationally constrained correlation matrix (actually the covariance matrix is used) in a multi-variate Gaussian distribution. Each random histogram is similarly integrated as in Eq. 1 to obtain a random copy of the α-age conversion curve. The dispersion of the 10 4 random samples of the conversion curve just gives the error bar at any given α value in Fig. 3. The computed results are given in Appendix Table 2. Table 1. Distance range used as criterion for member stars of each Gould's Belt protocluster and the numbers of removed foreground (n fg ) and background (n bg ) contaminating stars. Table 2. Data befind the plots: Columns 1-3 are for the average α histogram (after spline interpolation and normalized to an area of unity) and error bars (α, p and σp) in the left half of Fig. 2; columns 4-6 are for the α-age conversion curve and error bars (α, age and σage) in Fig. 3. The former curve is sampled at the bin centers of the histogram while the latter is at the borders. 2.1E-03 5.0 0.0E+00 0.0E+00 a The α = −2.25 histogram bin is an artificial one added to avoid the histogram from going to zero at α = −2. b The small positive histogram value of 6.2E-04 is due to the oscillation of the spline fitting function, which has little impact to this work.