20XX Vol. X No. XX, 000–000

11institutetext: Yunnan Observatories, Chinese Academy of Sciences, 396
Yangfangwang, Guandu District, Kunming, 650216, Yunnan, P. R. China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Chinese Academy of Sciences South America Center for
Astronomy, National Astronomical Observatories, Chinese Academy of Sciences,
Beijing, 100101, P. R. China
33institutetext: University of Chinese Academy of Sciences, Yuquan Road 19,
Shijingshan Block, Beijing, 100049, P. R. China
44institutetext: Departamento de Astronomía, Universidad de Chile, Las Condes,
7591245 Santiago, Chile
55institutetext: Instituto de Física y Astronomía, Universidad de Valparaíso,
ave. Gran Bretaña, 1111, Casilla 5030, Valparaíso, Chile
66institutetext: Centre for Astrophysics Research, University of
Hertfordshire, Hatfield AL10 9AB, UK
77institutetext: Xinjiang Astronomical Observatory, Chinese Academy of
Sciences, Urumqi 830011, China
88institutetext: Xinjiang Key Laboratory of Radio Astrophysics, 150
Science1-Street, Urumqi 830011, China
99institutetext: Shanghai Key Lab for Astrophysics, Shanghai Normal
University, 100 Guilin Road, Shanghai 200234, Peopleʼs Republic of China
Received 20XX Month Day; accepted 20XX Month Day

# Can mid-infrared spectral energy distribution quantitatively trace
protoplanetary disk evolution?

Mingchao Liu (刘明超) 112233 Jinhua He (何金华) 112244 Zhen Guo 5566 Jixing Ge (葛继兴)
7788 Yuping Tang (唐雨平) 99

###### Abstract

Infrared (IR) spectral energy distribution (SED) is the major tracer of
protoplanetary disks. It was recently proposed to use the mid-IR (MIR) SED
slope $\alpha$ defined between 2-24$\mu$m as a potential quantitative tracer
of disk age. We critically examine the viability of this idea and confront it
with additional statistics of IR luminosities and SED shapes. We point out
that, because the statistical properties of most of the complicated physical
factors involved in disk evolution are still poorly understood in a
quantitative sense, the only viable way is to assume them to be random so that
an idealized ‘average disk’ can be defined, which allows the $\alpha$
histogram to trace its age. We confirm that the statistics of the zeroth order
(luminosity), first order (slope $\alpha$) and second order characteristics
(concavity) of the observed MIR SED indeed carry useful information upon the
evolutionary processes of the ‘average disk’ and provide useful constraints to
future disk population synthesis modeling. We also demonstrate that intrinsic
diversities in MIR SED shapes and luminosities are always large at the level
of individual stars so that the application of the evolutionary path of the
‘average disk’ to individual stars must be done with care. The data of all
curves in plots are provided on GitHub111Disk-age package
https://github.com/starage/disk-age/.

###### keywords:

stars: pre-main-sequence — stars: protostars — infrared: stars — stars:
formation — protoplanetary disks

##  1 Introduction

The formation of low and intermediate mass stars has been intensively studied
and a standard picture of star formation (SF) has been there since a long time
ago (e.g., Shu et al. 1987; Larson 2003; McKee & Ostriker 2007; Evans et al.
2009). However, despite of tremendous efforts on observational studies (e.g.,
Myers et al. 1987; Palla & Stahler 2000; Haisch et al. 2001; Evans et al.
2009; Mamajek 2009, among many others), our quantitative knowledge of the star
formation process is still limited (e.g., Krumholz 2011). More quantitative
observational constraints on the earliest disk evolutionary stages are highly
desired to constrain disk evolution models, either at individual star level
(e.g., Zhu et al. 2010; Kimura et al. 2016; Kunitomo et al. 2020, 2021, among
many other models) or at population level (e.g., Emsenhuber et al. 2023), and
to test models of planet formation (e.g., Mordasini et al. 2009; Pignatale et
al. 2018; Burn et al. 2021; Raymond & Morbidelli 2022).

Although there have been many works to quantify disk evolution through the
observation of their host stars (see e.g., the review of Manara et al. 2023),
the most straightforward method should be using the IR emission of the disks
themselves. The first way is to utilise the IR colors in the color-color
diagram or color-magnitude diagram which is similar to the Hertzsprung-Russell
diagram (HRD; e.g., Palla & Stahler 2000) for stars. However, the colors were
usually used to identify YSOs, but not to quantify their evolution (e.g.,
Luhman et al. 2008; Rebull et al. 2010, among many others). The color-
magnitude diagram is also known to suffer from large uncertainty during the
earliest millions of years (Myrs) when the primordial disks exist (Soderblom
et al. 2014).

The second way is to use the characteristics of the IR SED. The most
frequently used parameter is the SED slope defined in logarithmic scale (known
as $\alpha$; Lada 1987; Adams et al. 1987; Greene et al. 1994). Since the
launch of the Spitzer Space Telescope (Werner et al. 2004), both the IR colors
and SED slope have been intensively used to explore the properties of nearby
populations of YSOs (e.g., Luhman et al. 2008; Rebull et al. 2010; Luhman et
al. 2010; Kirk et al. 2009; Megeath et al. 2012; Alves de Oliveira et al.
2012; Hsieh & Lai 2013; Broekhoven-Fiene et al. 2014; Dunham et al. 2015;
Azimlu et al. 2015, etc.). However, most of the existing works follow the
classical framework of disk evolution described in the IR SED classification
scheme of 0, i, f, ii and iii (Lada 1987; Adams et al. 1987; André et al.
1993; Greene et al. 1994). A typical example is to refine the lifetime
estimation of each SED classes using improved star counts (e.g., Wilking et
al. 1989; Evans et al. 2009).

In a recent work, Liu et al. (2023) proposed an innovative idea that the MIR
SED slope $\alpha$ defined between 2-24$\mu$m can be directly used to quantify
the disk evolutionary age in a novel manner not confined by the classical SED
classification scheme. Their key point is the striking similarity of the
$\alpha$ histograms among the 13 largest nearby protoclusters in Gould’s Belt
(GB). This also triggered them to propose that the large sample of disk SEDs
carry the information of disk evolution timescales and one only has the
opportunity to examine the star formation histories of individual
protoclusters after the common disk evolution timescales have been taken into
account. Indeed, only few works in literature paid attention to the peaks of
the $\alpha$ histograms, speculated that the peaks might reflect episodic SF
histories (e.g., Hsieh & Lai 2013). However, the similarity of the $\alpha$
histogram among all the GB protoclusters cannot be easily reconciled by
episodic SF events, because there is no known mechanism to synchronize these
parsec-scale events across so large a distance up to a kilo-parsec, which is
spanned by the entire Gould’s Belt.

With the belief that their average $\alpha$ histogram does trace disk
evolution, Liu et al. (2023) went further to divide the intrinsic disk
evolution processes into five stages, A, B, C, D and E, according to the
characteristics of the histogram. They also discussed how their staging scheme
can be better matched to known physical processes of disk evolution than the
traditional SED classes do.

From the theoretical point of view, it is reasonable to imagine an
evolutionary sequence of disk structure for a single YSO and it is fair to
allow the disk to have its own evolutionary age. As the major tracer of disks,
the IR SED or more specifically its SED slope $\alpha$, does have the
potential to serve as an independent age tracer for disks, just as the HRD
does for the central stars. However, as we will discuss below, the disk
evolution involves complicated physical processes and randomness. It is yet
uncertain how similar the SED evolution paths would be among YSOs formed in
different initial conditions and different local environments. The proposal of
Liu et al. (2023) is equivalent to assuming that these complication and
randomness can be effectively suppressed by the use of large star samples so
that their averaged $\alpha$ histogram mainly reflects the intrinsic disk
evolution. But a scrutiny on this point is necessary before we know to what
extent this disk-age tracer can be trusted. For this purpose, we first give a
brief review of the complicated disk evolution processes below.

###  1.1 The complicated relationship between protoplanetary disk evolution
and infrared SED

The relationship between the evolution of protoplanetary disks around low mass
stars and their MIR SEDs is expected to be complicated. In the earliest star-
formation stage, the accretion disk is deeply embedded in a gas envelope which
is the inner part of its natal molecular cloud dense core. The SED of such an
object is very red and usually too faint to be observable in IR due to high
extinction. The first opportunity to observe the youngest disk in IR is when
the stellar accretion has generated high speed bipolar jets to carve out a
pair of bipolar cavities in the envelope. The strength and shape of the
observed MIR SED in this stage sensitively depend on the inclination angle of
our sight line with respect to the bipolar cavities and the detailed
properties of the latter, and both thermal emission and scattered stellar
light can contribute (e.g., Kenyon et al. 1993; Whitney & Hartmann 1993;
Whitney et al. 2003). After then, due to the action of photoevaporation and
shock pressure of the jets, the bipolar cavities may expand quickly and the
envelope around the disk will be cleared (see some evidences in Jørgensen et
al. 2009; Furlan et al. 2016). The MIR SED in this stage still sensitively
depends on the inclination angle and evolves with time significantly.

Later on, the dust and gas in the disk may settle down into a layered
structure (e.g., Dullemond & Dominik 2004a). Cold dust grains in the middle
plane of the disk may aggregate into bigger ones to seed the formation of
protoplanets (e.g., D’Alessio et al. 2001, 2006), while the disk gas may
continue escaping from the surface layers through photoevapolation by FUV, EUV
and X-ray light from the central accreting star (Clarke et al. 2001; Alexander
et al. 2006; Rosotti et al. 2013) and maybe also due to magnetohydrodynamic
disk winds (e.g., Bai 2016). The inner disk edge may be puffed up by UV light
heating to form an inner wall to protect part of the disk against
photoevaporation through shadowing (Dullemond & Dominik 2004a). The MIR SED of
such a disk is determined by the thermal structure of the disk that is
regulated by various disk parameters such as the height, shape and thickness
of the inner wall and the flaring angle of the disk. The inclination angle may
still play some role. Formation of planets, particularly giant gaseous
planets, in this stage may induce dynamical features such as spirals, multiple
gaps of diverse widths or even a large inner hole in the disk (a transitional
disk; e.g., Alexander et al. 2006; Varnière et al. 2006; Gorti & Hollenbach
2009; Dodson-Robinson & Salyk 2011), which may have some impact to the MIR
SED. If a multiple star system is formed, the powerful perturbation of the
secondary star may truncate a significant part of the disk or even disperse
the entire disk earlier, resulting in a strong modification to the evolution
path of MIR SED (Kraus et al. 2012). Finally, the primordial disk will be
dispersed entirely by phototoevaporation so that the contribution of disk
emission to MIR SED diminishes.

In addition to the complexity, the disk evolution process also involves huge
randomness. For example, the initial envelope and disk masses, the initial
inner disk radius (depending on stellar mass), the extinction of parent
molecular cloud, the duration of accretion phase, and the number, masses and
radii of planets formed in the disk (see e.g., Robitaille et al. 2006) are all
random to some extent. The disk accretion to central star is also intermittent
(Kenyon et al. 1990), resulting in temporal variation of intrinsic MIR SED
shape and also additional uncertainty in observed SED shape when the different
photometric bands would not be observed simultaneously. Particularly, the
random inclination angle of the disk are believed to be able to alter observed
MIR SEDs significantly. (Whitney et al. 2003; Crapsi et al. 2008).

In this paper, we will first discuss how the above complexity would impact the
capability of the $\alpha$ histogram in quantifying disk evolution in Sect. 2.
Then, we report our scrutiny of the $\alpha$ histogram of Liu et al. (2023) in
Sect. 3. To strengthen our query into the reliability of the SED slope (a
first order characteristics of the SED) as a disk age tracer, we will further
employ the statistics on the second order SED characteristics (the IR SED
shape) and zeroth order SED characteristics (IR luminosities in the Spitzer
IRAC and MIPS $24\mu$m bands) of the same YSO samples in Sect. 4 and 5,
respectively. Particularly, we will be able to discern whether the most
important factor, the contamination from edge-on opaque disks, could prevent
the IR SED from being an effective age indicator. Finally, after a dedicated
discussion of the uncertainties in the $\alpha\sim-1$ histogram peak in Sect.
6, we draw our conclusions in Sect. 7.

##  2 Impact of uncertainties to the average SED-slope histogram

For many complicated processes and random factors discussed in Sect. 1.1, we
still do not have enough knowledge of them at the level of individual YSOs so
that their impact to MIR SED have to be considered as random processes at
population level. This is the case for the distribution of initial conditions
(e.g., initial envelope and disk masses, inner disk radius), total length and
intermittency of accretion phase, diversity of disk structures such as inner
wall, disk layers and flaring factor, the dust grain settling and growth,
planet formation and planet-disk interaction, extinction and observational
incompleteness. When a large population of YSOs are considered, as we are
discussing here, we can expect these random factors to be greatly suppressed
through averaging. However, concerning the average $\alpha$ histogram of Liu
et al. (2023, their Fig. 1; also see an improved version in our Fig. 4 in
Appendix), the sample sizes of the earliest stage A, the earlier part of Stage
B and the latest stage E are not large enough so that the advantage of
averaging does not work. Fortunately, the earliest stages are brief so that
they only have limited impact to age estimation in later stages. Furthermore,
we should also remain alert in future on possible distortion of the $\alpha$
histogram by any yet unrecognized non-negligible parameter correlations.
Below, we briefly discuss several better understood factors among them.

The intermittent disk accretion has been studied since a long time ago (e.g.,
Kenyon et al. 1990; Hartmann & Kenyon 1996; Megeath et al. 2012; Hartmann et
al. 2016) and has become a hot topic (e.g., Vorobyov & Basu 2015; Caratti o
Garatti et al. 2017; Contreras Peña et al. 2017; Guo et al. 2020, 2024a,
2024b; Park et al. 2021). But the temporal behavior of such bursts is still
not well predictable (Fischer et al. 2023). There have been clues indicating
that the intermittency could be related to random gas in-fall onto the disk
(Kuffmeier et al. 2018). Therefore, currently, the best choice is to consider
the accretion bursts as a random process.

Another important complication, the gaps, central hole or other dynamical
effects created in disks by protoplanets or stellar companions, can also be
treated as random. Despite of a long history of intensive modeling and
observations of planet formation and its interaction with the disk (e.g.,
Pollack et al. 1996; Alibert et al. 2005; Rosotti et al. 2013; Baruteau et al.
2014; Dong et al. 2015a, b; Andrews et al. 2018; Long et al. 2018; Emsenhuber
et al. 2021, among many others), the number of planets, their masses and orbit
radii, or similar paramters of stellar companions, are still hard to predict
for a real population of protoplanetary disks.

Differential extinction can modify observed MIR SEDs. Generally, only
foreground extinction of interstellar medium can be easily corrected in
observations. The remaining, more localized extinction can be classified into
two categories: 1) the one due to local cloud and circumstellar envelope and
2) the one due to disk. The former is known to be larger in denser clouds
where the younger YSOs of Stages A and B reside (e.g., Getman et al. 2014;
Azimlu et al. 2015, among many others). However, in any evolutionary stage of
disk or in any given cloud region, its value for different YSOs is more or
less random. Beside the randomness, the remaining deterministic effect of MIR
SED reddening indeed has only limited impact to the $\alpha$ histogram due to
the paucity of objects in these early stages.

The second category of extinction (by disk) sensitively depends on the
geometrical and optical thickness of the disk and some secondary parameters
such as flaring and self-shadowing by inner wall. Its impact is expected to be
large among Class 0, i and f YSOs in the earliest Stages A and B when the disk
is geometrically thick. Even for later Class ii objects in Stages C and D,
Crapsi et al. (2008) had suggested that the impact of disk inclination angle
would be as high as $30\%$. However, we would like to figure out an intriguing
difference between the inclination angle effects and the contamination from
edge-on opaque disks (opaque disks that totally block out the inner hot disk
region): while the former is likely always important, the latter is not
necessarily so. For example, if the disks are always very opaque during most
of their lifetime, as demonstrated by the model grid of Robitaille et al.
(2006), edge-on disks may become so faint at near IR (NIR), due to strong
extinction, that they almost never appear in our Spitzer/IRAC selected sample.
In this case, the edge-on opaque disks do not contribute to modifying observed
$\alpha$ values, while the inclination angle effect is still important in
determining the fraction of observable disks (depending on the cavity opening
angle) and in reddening the observed SEDs to a smaller extent (depending on
structural details of the disk surface) in all evolutionary stages. As we will
further explore in more detail in the next sections, the real impact of edge-
on disk contamination is indeed not dominant in the GB samples.

Observational incompleteness of YSOs at population level can have strong
impact to $\alpha$ histograms too. The main sources of incompleteness are the
confusion with background/foreground main sequence stars (e.g., Harvey et al.
2007) and high extinction. The former mainly affects the Class iii objects in
Stage E due to the similarity of their SEDs to that of stellar photosphere,
while the latter leads to the missing of YSOs with lower luminosity. This fact
reinforces the vulnerability of the average $\alpha$ histogram in the earliest
and latest stages A, early B, and E. Perhaps some on-going projects such as
SESNA (Pokhrel et al. 2020) could improve it to some extent with a better way
of source extraction in the Spitzer maps.

In a summary, while the randomness in the disk evolution processes could be
suppressed by using large samples of YSOs, the relationships between disk
parameters are still poorly known and there is some risk to treat all the disk
parameters as independent random factors. Thus, further testing of the
reliability of the SED slope as a disk age tracer will be beneficial.

##  3 Scrutinize the MIR SED-slope histogram

As the first step of our scrutiny, we closely examined the methodology of Liu
et al. (2023). Although their judgement of cluster membership, new extinction
correction, handling of Spitzer MIPS $24\,\mu$m saturation cases and the
definition of MIR SED slope $\alpha$ are acceptable, we identify two important
improvements. The details of the improvements, an updated version of the
average $\alpha$ histogram and a simple discussion of its potential to trace
disk evolutionary age are given in Appendix A. The main conclusion is that it
is mainly the Stages B, C and D of Liu et al. (2023) that potentially have the
power to trace disk ages.

We also compare the observed $\alpha$ histogram to that of the disk model grid
of Robitaille et al. (2006) in Appendix A.2 in order to verify whether the
former contains useful information of disk evolution to constrain disk models.
Despite of a number of important limitations of the model grid (Sect. 2.2.4 of
Robitaille et al. 2006; Robitaille 2008), it has integrated the disk physics
known at that time. It is found that the observed and the modeled $\alpha$
histogram are very different. The observation $\alpha$ histogram can provide
useful constraints to disk model grids.

##  4 Test the evolution of MIR SED shape and its diversity

Real SEDs usually take curvy shapes, which is jointly determined by the
multiple factors discussed above. If the SED slope $\alpha$ really can trace
disk evolution, we should expect that the change of the average SED shape
across the bins of the $\alpha$ histogram mainly reflects the evolution of
disk structure from a spherical envelope to bipolar cavities to a
geometrically thick disk and finally to a geometrically thin disk, while the
diversity of the SED shapes measures the randomness of the physical factors.
If the contamination from edge-on opaque disks would have played a major role
during any of the Stages B, C or D, we should expect more concave SED shapes
and an enhancement of the SED-shape diversity at larger $\alpha$ value, due to
higher extinction at shorter wavelengths and due to the contribution of
scattered light in optical and NIR, according to the simulations of Whitney et
al. (2003); Robitaille et al. (2006) etc. We can verify these trends among the
same sample of GB YSOs.

To character the SED shape, we extend the linear fit to the SEDs when defining
the slope $\alpha$ to a quadratic one

$lg(\lambda F_{\lambda})=alg^{2}(\lambda)+blg(\lambda)+c$ (1)

so that we can use the second order coefficient $a$ to characterize the shape
(or concavity) of an SED. An SED is overall convex when $a<0$, lacking second
order concavity when $a=0$, overall concave when $a>0$. An example of each
case is shown in Fig. 1. Note that this concavity does not reflect the
Silicates feature around $9.8\mu$m. The coefficient $b$ is similar to
$\alpha$, but we still use the latter to represent the SED slope, because it
is more consistent with the tradition in literature.

The SED shape reflects the joint effects of the wax and wane of the major
emission components (e.g., the central star, inner and outer parts of the
disk, bipolar cavities and envelope) and extinction. Previous works used SED
slopes or colors over diverse wavelength ranges such as $3.6-8.0\mu$m and
$5.8-24\mu$m (e.g., Fang et al. 2009, 2013) or a turnoff wavelength
$\lambda_{\rm turnoff}$ of the SED (Cieza et al. 2007) to characterize the SED
shape. The second order coefficient $a$ in our Eq. 1 can quantify such
information in the simplest manner enabling statistical study.

Figure 1: Examples of quadratic polynomial fitting (full line) of observed
SEDs (dots) to illustrate how the second order coefficient $a$ quantifies the
SED concavity. The observed SEDs have been shifted to the same value of zero
in Ks band (left most wavelength) for clarity. Figure 2: The average
($\overline{a}$, black, the error bar is the uncertainty) and sample standard
deviation ($\sigma_{a}$, gray) of the SED-shape parameter $a$, as defined in
Eq. 1, for each $\alpha$ bin of width $\Delta\alpha=0.2$ (left for all the
5194 observed samples in Gould’s Belt; right for the model grids of Robitaille
et al. 2006, with SED averaged over the 50 model apertures). Only the bins
with observed star number $\geq 10$ are shown. The long vertical gray dashed
lines mark the border lines of the disk evolution Stages B, C, D and E of Liu
et al. (2023). The horizontal gray dashed line at $Y=0$ is for comparison. The
data of all the curves are given on the GitHub mentioned in Abstract.

The average ($\overline{a}$) and standard deviation (STD; $\sigma_{a}$) of the
parameter $a$ in each $\alpha$ bin of width $\Delta\alpha=0.2$ are shown
together in the left panel of Fig. 2. The former quantifies the average second
order shape of SED while the latter characterizes the SED diversity. The error
bar of $\overline{a}$ is simply $\sigma_{a}$ divided by the square root of
star number. The most striking feature is that $\sigma_{a}$ nearly
monotonically decreases with decreasing $\alpha$, which quantitatively
reproduces the same overall trend of diminishing SED diversity revealed in
Fig. 6 of Liu et al. (2023). This trend uncover the fact that the SED
diversity (perhaps also disk-structure diversity) continues to diminish with
disk evolution. Furthermore, the nearly linear overall correlation between
$\sigma_{a}$ and $\alpha$ hints that the dominant causes of the SED diversity
must have some simple relationship with disk thermal structure (traced by
$\alpha$) in the entire lifetime of a protoplanetary disk. Further detailed
modeling such as those of Kimura et al. (2016); Schib et al. (2021); Karam &
Sills (2023) is needed for a quantitative interpretation. Distinct from the
STD, $\overline{a}$ shows a broad peak in the range of $\alpha<0$ but remains
around zero for $\alpha>0$. We will interpret these features stage by stage
below. These features indicate that the MIR SEDs are not dictated purely by
random factors, but by the physical relationships of the disk parameters.
However, another valuable fact is that $\sigma_{a}\geq\overline{a}$ is roughly
true for almost all $\alpha$ values. Even around the $\alpha\sim-1$ histogram
peak, $\sigma_{a}$ is still comparable to the average. This means the SED
diversity is always significant in all disk evolutionary stages, which
cautions us that the average view of disk evolution in terms of MIR SED may
not be simply applicable to individual disks.

Stage B is characterized by very high SED diversity ($\sigma_{a}$) but nearly
zero average SED concavity ($\overline{a}$). This fact simply expresses that
there does not exist representative SED in this stage and half of the SEDs are
convex and half concave. Because the YSOs in this early stage very likely
still possess a thick envelope with bipolar cavities, inclination angle should
be the most sensitive parameter to determine the observed SED shape. However,
as we mentioned in the beginning of this subsection, if the contamination from
edge-on disks is the dominant factor controlling the observed $\alpha$
histogram, we should expect a large fraction of YSOs to be edge-on opaque
disks that have moved here from Stage C (because of many more stars in Stage C
than in B) and they have more concave SED shapes, according to simulations
such as Whitney et al. (2003, see e.g., their Fig. 4). Thus, in Fig 2, the
average SED concavity $\overline{a}$ is expected to be larger in this stage
than in C, which is not the case in the figure. Therefore, we can conclude
that the member stars in this stage should not be dominated by the
contamination from edge-on opaque disks modeled by Whitney et al. (2003).
Instead, it can be interpreted by the gradually clearing of the spherical
envelope. We also notice that the $\overline{a}$ curve suffers from large
oscillations when $\alpha>1$ but converges better to zero in the range of
$1>\alpha>0$. This can be understood as the consequence of the increase of
star number with decreasing $\alpha$, as can be seen in the upper panel of
Fig. 4. This interpretation is also consistent with the larger error bars in
those oscillating $\alpha$ bins in Fig. 2.

Stage C is the right half of the sharp $\alpha\sim-1$ peak. While $\sigma_{a}$
continues to linearly drop with decreasing $\alpha$ in this stage, the average
SED concavity $\overline{a}$ first rises up quickly to a peak value of $\sim
0.95$ around $\alpha=-0.7$ (i.e., the SEDs become more concave) and then falls
down slowly (i.e., become slightly less concave) towards the left border at
$\alpha=-1$. Similarly, the $\alpha>-0.7$ part of this stage cannot be
dominated by the contamination from edge-on opaque disks moved here from the
$\alpha<-0.7$ part, because otherwise we should see more concave SEDs (with
larger $\overline{a}$) in the former part than in the latter, which is not the
case in the figure. The sharp rise of $\overline{a}$ with decreasing $\alpha$
can be interpreted by physical reasons – very likely the stabilization of the
disk to a structure maintaining significant $24\mu$m excessive emission but
increasingly less at NIR. This is possible if a large part of the inner disk
is gradually shadowed by the building up of a puffed-up inner rim (Dullemond
et al. 2001; Dullemond & Dominik 2004b, a) while the outer disk is still
flaring (remaining warm by illumination to maintain enhanced $24\mu$m
emission). The establishment of the disk shadow drives the slow decrease of
$\alpha$, while the stabilization of disk structure means the slow down of
disk evolution, which explains the sharp rise of the $\alpha$ histogram
towards the peak at $\alpha\sim-1$. Then, given that the $\alpha>-0.7$ part of
this stage is already not dominated by the contamination from edge-on opaque
disks, as discussed above, the $\alpha<-0.7$ part that has many more stars
should be even less so. The slow decrease of the average SED concavity
$\overline{a}$ with decreasing $\alpha$ in this part also can be interpreted
by physical processes – possibly the slow photoevaporation of the flaring
outer disk and/or gap opening by proto-giant-planets, both of which have the
capability to reduce the $24\mu$m emission so as to decrease both the SED
slope $\alpha$ and the concavity $a$ to some extent.

Stage D is the left half of the $\alpha\sim-1$ peak. Both the average
$\overline{a}$ and dispersion $\sigma_{a}$ of the SED concavity have stopped
the decreasing trends with decreasing $\alpha$. There is even the indication
of a small bump in both curves in the range of $-2<\alpha<-1.4$, which means
the MIR SEDs become slightly more concave and diverse than expected from the
trends in Stages B and C. It is obvious that this stage is not dominated by
edge-on disk contamination, given that none of Stages B and C is. The bump can
be interpreted also by physical reasons such as the creation of a large inner
hole or wide gaps by giant planets in some disks (the transitional disks) or
any yet unknown special properties of the final disk dispersal processes (very
likely by photoevaporation). In addition, the contamination of the sample by
some debris disks is also possible.

In Stage E, both $\overline{a}$ and $\sigma_{a}$ decrease with decreasing
$\alpha$ and $\sigma_{a}>\overline{a}$ is always true. These trends should be
discuss in future when the sample would be more complete in this stage.

We also introduce the disk-SED model grid of Robitaille et al. (2006) to
compare with the observed SED concavities, as we did in the comparison with
the observed $\alpha$ histogram in Sect. A.2. We fit the SED concavity $a$ of
all their disk models and show the average and dispersion of the model SED
concavities as functions of MIR SED slope $\alpha$ in the right panel of Fig.
2. It can be seen that the distributions of both the average and dispersion of
the model concavities are very different from the observed ones in the left
panel, which means the of the model disk parameter distributions can be
constrained to some extent by comparing them with the observed ones.

In a summary, the MIR SED concavity $a$ has allowed us to make a conservative
point that, although the inclination angle effect could still be important in
all of the Stages B, C and D (because of the always large SED diversity), the
contamination from edge-on disks is never the dominating factor in any of
these stages. This suggests that most of those edge-on disks are likely opaque
even at IR wavelengths so that most of them actually did not enter the
Spitzer/IRAC selected sample. Both the MIR SED slope $\alpha$ and concavity
$a$ are consistent with the intrinsic evolution of the disk structures in
these stages, including the clearing of the spherical envelope in Stage B, the
formation and photoevaporation of shadowed flaring disks in Stage C and the
development of transitional disks and the final runaway disk dispersal in
Stage D.

##  5 Test the monochromatic luminosity trends

Figure 3: The distribution of average luminosity $\log_{10}(\nu F_{\nu}D^{2})$
vs MIR SED slope $\alpha$ (in each $\alpha$ bin of width 0.2) in the 3.6$\mu$m
(left) and 24$\mu$m (right) bands, respectively, for the Gould’s Belt YSOs
(red) and the 200710 disk models of Robitaille et al. (2006, green; averaged
over the 50 model apertures), where $D$ is the average GAIA distance to each
observed protocluster and 250 pc for the models. The small error bars show the
uncertainty of the average luminosity while the large error bars (cyan for
observed YSOs and orange for the models) are the sample dispersion of the
luminosities.

Our next attempt is to examine the statistics of the zeroth order SED feature,
the monochromatic IR luminosities of the YSOs. We will discuss whether the
important factor, the edge-on disk contamination, is a dominating factor and
whether the results are consistent with disk structure evolution. According to
the modeling of Wood et al. (2002); Whitney et al. (2003), edge-on opaque
disks are not only redder, but much fainter at NIR wavelengths due to the high
extinction of the optically thick disk. If any sub-region of the $\alpha$
histogram contains a significant fraction of edge-on disks, they may stand out
as objects with significantly lower IR luminosity than the others and the
average luminosity in that sub-region will be lowered.

For this purpose, we plot the luminosities of the GB YSOs in the IRAC and MIPS
$24\mu$m bands as functions of the SED slope $\alpha$ in the left column of
the Fig. 5 in Appendix B. The luminosity is represented by $\nu F_{\nu}D^{2}$
where $\nu$ is the frequency, $F_{\nu}$ is the flux and $D$ is the averaged
GAIA distance (weighted by signal-to-noise ratio) to each protocluster derived
similarly as in Liu et al. (2023). YSOs without reliable $24\mu$m flux are
excluded because the version of $\alpha$ used in this work cannot be well
defined without it. To facilitate discussion, we also compute the average of
$\log_{10}(\nu F_{\nu}D^{2})$ in each $\alpha$ bin of width 0.2 (red curves)
in the Stages B, C and D which we are interested in. Two sets of error bars
are shown, with the larger one (cyan) being the sample dispersion and the
smaller one (red) the uncertainties of the averages (dispersion divided by the
square root of star number; both defined also in the logarithmic space).

In order to better understand the statistical behaviors of the observed YSOs,
we also repeat the above operations to the model SEDs of the disk model grid
of Robitaille et al. (2006). We compare the average and dispersion of
luminosities between the models and observations in the second column of Fig.
5.

To facilitate discussions here, we only pick out the most interesting
$3.6\mu$m and $24\mu$m luminosity curves and show them side by side in Fig. 3.
Interesting trends can be seen in the two plots. The observed $3.6\mu$m
luminosities (red curve) show an overall increasing trend with decreasing
$\alpha$. However, the trend becomes flat in Stages C and D. Oppositely, the
$24\mu$m luminosities show an overall decreasing trend with decreasing
$\alpha$. The decreasing trend in Stages C and D is slightly steeper than in
Stage B. Another interesting point is that the larger sets of error bars are
almost homogeneous in all $\alpha$ bins despite of the huge difference in star
numbers among the bins, which hints that the large scatter in the plots are
mainly not due to random measurement errors, but due to intrinsic physical
diversities of the YSOs in all three stages. Possible contributors to the
luminosity diversity include different stellar and disk masses and
intermittent stellar accretion. Our analysis of the apparent fluxes (not
plotted here) shows that the detected YSOs all have their IRAC and MIPS
$24\mu$m fluxes far above the nominal sensitivity limits of the Spitzer bands
given by Werner et al. (2004), which means that the luminosities in Fig. 3 do
not suffer from sensitivity limitation.

These results do not support the contamination from edge-on disks to be
important in the sample. At the level of individual stars, we do not see the
expected subsample of YSOs standing out as distinctive low luminosity sources
in any of the three stages. Although we cannot exclude the existence of some
edge-on opaque disks, they must not be numerous in the observed sample. At the
population level, if Stage B is strongly contaminated by edge-on opaque disks
moved there from Stage C where there are many more member stars than in B, due
to the strong disk extinction, we should expect the Stage B YSOs to have
significantly lower luminosities than those in C. However, at least the
$24\mu$m luminosities in the right panel of the figure show an opposite trend
which should be interpreted by stronger excessive $24\mu$m emission from the
more massive cool spherical envelope in Stage B. The grand anti-correlation
between the $3.6\mu$m luminosity and $\alpha$ in Stage B can be explained by
either the dominance of edge-on opaque disk contamination or the evolution of
the YSOs (i.e., progressively clearing the envelope) or both, because both
factors function through the change of opacity (reddening). However, only the
latter is consistent with the trend in the $24\mu$m luminosities. Similarly,
were Stage C dominated by edge-on disks moved there from Stage D due to strong
extinction of edge-on disks, the YSOs in this stage should have significantly
lower luminosities than those in Stage D, but this is not the case in either
panel of the figure.

After excluding the contamination from edge-on disks, we shall examine whether
the luminosity trends in Fig. 3 are consistent with disk evolutionary
processes. Both the $3.6\mu$m and $24\mu$m luminosity trends in Stage B have
been reasonably suggested in the previous paragraphs to be driven by the
clearing of the envelope. Further on, the comparable $3.6\mu$m luminosities in
Stages C and D can be due to the slow evolution of the stellar luminosity,
because stellar emission has become the major contributor in these stages. The
more steeply decreasing $24\mu$m luminosities with decreasing $\alpha$ in
these two stages can be explained by the formation of giant planets that open
broad gaps to destroy the part of disks that is mainly responsible for the
$24\mu$m emission or by the photoevaporation of the outer disks. It is worth
noting that the transition point at $\alpha=-0.3$ in both panels of the figure
is close to the $\alpha$ value where the average MIR SED concavity
$\overline{a}$ rises up quickly from 0 to $\sim 1$ as the $\alpha$ decreases
in Fig. 2, which implies that they may be two manifestations of the same
process: the start of giant planet formation.

The luminosity curves of the disk model grid of Robitaille et al. (2006) in
Fig. 3 (green curves) show increasing trends with increasing $\alpha$ value
and a peak around $\alpha=0.4$ in both 3.6 and 24$\mu$m bands. The trends are
in opposite sense than the observed 3.6 $\mu$m luminosity distribution and
significantly steeper than the observed 24 $\mu$m luminosity trend. In
addition, the sample dispersion of the model grid (error bars in orange) are
definitely larger in the $\alpha>0$ range than in the $\alpha<0$ range, which
is also distinct from the uniform sample dispersion (error bars in cyan) of
the observed luminosities. The large differences indicate that the
luminosities can provide meaningful constraints to disk model grids as well.

In a summary, the monochromatic IR luminosities do not support the
contamination from edge-on opaque disks to be a dominant factor in Stages B, C
and D in the $\alpha$ histogram. Instead, they are consistent with disk
evolution. This conclusion casts some doubt on the high fraction of more than
$30\%$ of mis-classified Class ii edge-on disks in the model of Crapsi et al.
(2008). Possible improvements to their model are to consider round inner disk
wall as suggested by Isella & Natta (2005), to consider optically thicker
disks so that edge-on disks may become too faint to enter the observed sample,
or alternatively to consider more transparent disks so that many edge-on disks
have little impact to MIR SEDs, or to assume a geometrically thinner disk. On
the other hand, the nearly uniform, large, luminosity scatter in all stages in
Fig. 3 also warn us that the $\alpha$ histogram is mainly suitable for the
‘average disk’ while applying it to individual stars must take into account
the large intrinsic diversity.

##  6 Large uncertainty in the SED-slope histogram peak?

The $\alpha$ histogram peak around $-1$ mimics a Gaussian function to which
one cannot help suspecting that it might be still dominated by random factors,
instead of by disk evolution. Particularly, the large intrinsic diversity of
the disks revealed in Figs. 3 and 2 hints that the sample dispersion may not
be thoroughly compressible through averaging. However, if we assume that the
peak shape is still dominated by randomness, the underlying true $\alpha$
histogram of the average disk would be even narrower than currently observed.
This hints at a more extreme disk evolution pattern in which the MIR SED slope
$\alpha$ of a given YSO must evolve from a larger value (e.g., 0) into the
peak value of $\sim-1$ in a sudden jump, which is not likely for huge objects
like the protoplanetary disks. Furthermore, the slight shift of the histogram
peak from $\alpha=-1$ to the asymmetric position $\alpha=-1.03$ in the
improved version of the $\alpha$ histogram in Fig. 4 also hints that this peak
is not a perfect Gaussian function.

##  7 Summary

We have critically examined the uncertain factors behind the histogram of MIR
SED slopes $\alpha$ of Gould’s Belt YSOs. We have to confess that still little
is known of the statistical nature of most of the complicated physical factors
involved in the evolution of proto-planetary disks so that they have to be
assumed to be random factors for the moment. It is obvious that the $\alpha$
histogram cannot trace disk evolution well in the earliest Stage A and the
latest Stage E due to sample incompleteness. By using additional information
from monochromatic luminosities and SED concavity, we are able to show that,
although the disk inclination angle might be always important, the
contamination from edge-on opaque disks is never a dominating factor to the
observed MIR SED at a population level. The statistics of the luminosity, SED
slope $\alpha$ and concavity can trace the overall evolutionary sequence of an
‘average disk’ that includes envelope clearing, formation of self-shadowed
flaring disk, gap or hole opening by planet formation, photoevaporation etc. A
linear decrease of MIR SED diversity with decreasing $\alpha$ is found (left
panel of Fig. 2), but needs further modeling to explain quantitatively.
However, the application of these results to individual stars should be done
only with an adequate account for the large intrinsic diversities of the IR
tracers that are not necessarily purely random in nature.

###### Acknowledgements.

ML jointly motivated this work, analyzed the data and drafted the manuscript;
JH jointly motivated this work, participated in the data analysis and revised
the manuscript; the other co-authors discussed and improved many parts of the
work and manuscript. This work is supported by the Natural Science Foundation
of Yunnan Province (No. 202201BC070003). ZG is supported by the ANID FONDECYT
Postdoctoral program No. 3220029. ZG acknowledges support by ANID, –
Millennium Science Initiative Program – NCN19_171. This work is sponsored (in
part) by the Chinese Academy of Sciences (CAS), through a grant to the CAS
South America Center for Astronomy (CASSACA) in Santiago, Chile.

##  Appendix A Improvements to the histogram of SED slopes

###  A.1 The improvements

We have carefully examined the procedures of Liu et al. (2023) and found that
the following tow aspects deserve improvements: 1) the resolution of the
average $\alpha$ histogram is far from optimal; 2) their choice of averaging
the $\alpha$ histograms among the 13 protoclusters is only one of the possible
schemes to derive the average $\alpha$ histogram, while the impact of
different schemes to the resulting average $\alpha$ histogram is still
unclear.

There are several empirical approaches to determine the optimal histogram bin
width of specific kinds of sample distributions. Usually the optimal bin width
can be smaller for a larger sample size. If the random samples are averaged in
sub-groups before final averaging, the effective sample size could be slightly
smaller than the total sample size, depending on the weighting among sub-
groups. However, this shrinkage of sample size is not significant and can be
neglected, according to our practice on the sample in this work. Therefore,
the choice of the ‘Freedman-Diaconis’ rule (Freedman & Diaconis 1981) which is
good enough for a singly peaked histogram can be justified for the work of Liu
et al. (2023). However, they arbitrarily adopted a conservative bin width of
$\Delta\alpha=0.5$ without any explanation. If we adopt the total size of 5194
YSOs in all the 28 GB protoclusters surveyed by Spitzer Space Telescope, we
derive the optimal bin width of $\sim 0.082$. As can be seen from the top
panel of Fig. 4, the higher $\alpha$ resolution still produces a smooth
histogram but presents more details of the distribution.

There are also different ways to make the $\alpha$ histogram for a large
sample of YSOs. The most straightforward way is to merge all the YSOs of the
28 GB protoclusters surveyed by Spitzer Space Telescope into a single sample.
A drawback of this approach is that large protoclusters, such as Orion A, may
dominate the resulting histogram so that it may be biased by the special
properties (e.g., the temporal fluctuation of star formation rate) of few
largest clusters. Liu et al. (2023) adopted the natural star grouping of the
individual GB protoclusters and averaged the histograms over the 13 largest
ones to improve the imbalance of cluster weights. However, the 13 clusters
still have a broad range of sizes from 70 to 2435 member stars. Some of the
large protoclusters can be broken into smaller sub-clusters. Therefore, we
visually assess the clustering of member stars in each of the 13 GB
protoclusters in the sky plane and also in their distribution of $\alpha$
value and further divide them into sub-clusters (36 in total; the details will
be given in a future paper). These sub-clusters still show similar $\alpha$
histograms, but with the differences among them being larger than among the 13
protoclusters.

We compare the resulting (averaged) $\alpha$ histograms from the three
approaches in the upper panel of Fig. 4 to assess the impact. The histogram
made by taking all the GB YSOs as a single sample (cyan) shows a peak position
at a slightly larger $\alpha$ value than the other two, very likely because it
is dominated/biased by the largest cluster, Orion A, which has contributed
nearly one half of member stars in the ‘single sample’ case. The histograms
from the 13 largest GB protoclusters (red) and from the 36 sub-clusters
(black) are more similar in overall shape, which hints that more deliberately
grouping the YSOs may not benefit further. However, there still are small but
potentially important differences between the latter two: 1) the 36-sub-
clusters histogram shows a stronger and narrower peak around $\alpha=-1.03$
which best informs an intriguing asymmetry in the peak shape and 2) shows a
slightly stronger secondary peak around $\alpha=-2.5$. Because the grouping of
the 36 sub-clusters has given the diverse sub-clusters the most homogeneous
weights in the averaging, we consider it as the best approach.

Similar as Liu et al. (2023), we fit a cubic spline function (black dashed,
upper panel of Fig. 4) to the best $\alpha$ histogram using the function
splrep of the Python package scipy.interpolate, with an empirical choice of
the smallest possible smoothing parameter $s=0.08$ to suppress unrealistic
fluctuations in the fitted curve. Then, we integrate the fitted $\alpha$
histogram to derive an $\alpha$-to-age conversion curve, as shown in the
bottom panel of Fig. 4 (black curve). Compared to the version of Liu et al.
(2023, red curve), our new curve is steeper around the $\alpha=-1$ peak and
extends down to the incomplete range of $\alpha=-3\sim-2$.

Figure 4:  The upper panel shows the $\alpha$ histograms derived with three
star grouping schemes: 1) taking all the 5194 GB YSOs as a single sample
(cyan); 2） averaging the histograms over the 13 largest GB protoclusters (red;
4909 YSOs in total); 3) averaging the histograms over the 36 sub-clusters
(black; 4312 YSOs in total). The black dashed curve is a smoothed cubic spline
fitting to the best histogram, the 36-sub-cluster case. The gray shaded
histogram is the $\alpha$ histogram of the 200710 disk models of Robitaille et
al. (2006, SED averaged over the 50 model apertures). The vertical dashed
lines show the five disk evolution stages suggested by Liu et al. (2023). The
bottom panel shows the new $\alpha$-to-disk age conversion curve (black)
derived from the best $\alpha$ histogram in the top panel, which can be
compared to the result from Liu et al. (2023, red curve). The vertical dashed
line mark the critical value of $\alpha=-2$ below which our sample becomes
incomplete. The data of the smoothed $\alpha$ histogram and the new alpha-to-
age conversion curve are given on the GitHub mentioned in Abstract.

###  A.2 The new SED-slope histogram

Now we discuss the new features of the updated $\alpha$ histogram and re-
examine its potential to trace disk age. First of all, if we adopt the same
criteria of Liu et al. (2023) to divide evolutionary stages according to the
characteristics of the $\alpha$ histogram, the five evolutionary stages, A, B,
C, D and E as proposed by them and shown in the upper panel of Fig. 4 still
roughly agree with the statistical features of the new $\alpha$ histogram.
Thus, we do not tend to update it.

The Stage A ($\alpha>2.5$) involves very few observed YSOs. The histogram in
this stage is not reliable and thus cannot trace the evolutionary age well.
Fortunately, however, this stage is likely very brief in nature so that its
impact to the estimation of disk age in more evolved stages is limited.

In Stage B ($2.5\geq\alpha>0$), our new $\alpha$ histogram shows a clearer
ramp-up with decreasing $\alpha$ than in that of Liu et al. (2023). The nearly
linear shallow distribution may allow the $\alpha$ histogram to trace disk
evolution in a simple linear manner. However, the total number of YSOs
involved in this stage is still small so that the impact of large
uncertainties still need to be further discussed (see in next section).

Stages C ($0\geq\alpha>-1$) and D ($-1\geq\alpha>-2$) bracket the major peak
of the $\alpha$ histogram. However, with the finer resolution, this peak is
definitely sharper and taller than in that of Liu et al. (2023). In addition,
we also notice a small shift of the peak position to $-1.03$, which indicates
some asymmetry in the shape of the peak. If confirmed by future observations,
the asymmetry would indicate that the clearing of the primordial disk is
faster than establishing it.

The peak of the $\alpha$ histogram around -1 has an interesting indication. It
is known that a perfectly flat disk heated by UV radiation from an extended UV
emission region around the central star or dominated by steady disk-accretion
heating has a theoretical MIR SED slope of $-1.33$ (Dullemond et al. 2007).
The observed $\alpha$ peak of $-1$ (or more exactly $-1.03$) is significantly
larger, which means the disk structure in the most stable evolutionary phases
must have non-negligible flaring degree. This conclusion also agrees to later
observations (e.g., Avenhaus et al. 2018). Dullemond et al. (2007) showed that
a disk model with an MIR SED slope $\alpha=-1$ should assume a radial
temperature profile of $T_{eff}\propto r^{-\frac{2}{3}}$, which is shallower
that of a flat disk model, $T_{eff}\propto r^{-\frac{3}{4}}$.

Concerning $\alpha$ as a disk evolution tracer, even if this is true, the
sharp peak of its histogram indicates that the disk must evolve very slowly
near the peak so that the $\alpha$ values within the full width at half
maximum (FWHM) of the peak ($-1.3<\alpha<-0.7$) might not be a sensitive
tracer of the disk evolution. Nevertheless, as we will show in Sect. 5, even
within the FWHM of the peak, $\alpha$ still has the potential to trace the
disk evolution to some extent.

The completeness of samples in Stage E is not great because only few
protoclusters have their older Class III member stars (Liu et al. 2023)
successfully identified, so it is difficult for $\alpha$ to accurately reflect
the age of disk evolution. Due to the same reason, we also do not discuss the
secondary peak around $\alpha=-2.5$.

To check whether the observed $\alpha$ histogram can be used to constrain
model grids in a statistical sense, we compare it to the disk model grid of
Robitaille et al. (2006). We compute the MIR SED slope $\alpha$ for all the
200710 disk models in the same manner as for the observed SEDs and show their
histogram in the upper panel of Fig. 4 (gray filled). Its striking difference
from the observed one demonstrates the usefulness of the observed $\alpha$
histogram to constrain future disk evolution model grids.

In a summary, solely judged from the properties of the derived average
$\alpha$ histogram itself, the MIR SED slope $\alpha$ may not be a reliable
tracer of disk evolution in the earliest Stage A, the earlier part of Stage B,
the last Stage E, and even in the core part of the $\alpha=-1.03$ peak.
However, Stages B, C and D of the histogram still can provide useful
constraints to future disk evolution model grids.

##  Appendix B All luminosity plots

The correlations between all the observed monochromatic luminosities in the
four IRAC and the MIPS 24$\mu$m bands with the corresponding MIR SED slope
$\alpha$ are presented in the left most column of Fig. 5. A comparison with
similarly computed monochromatic luminosities of the disk-SED model grids of
Robitaille et al. (2006) is given in the second column of the figure. See the
detailed discussions of them in Sect. 5. For completeness, we also plot
similar correlations between the luminosities and the SED concavity $a$ among
the observed YSOs and among the models of Robitaille et al. (2006) in the
third and fourth columns of the Fig. 5. They may be useful for constraining
future model grids. We do not discuss the details here.

Figure 5: The distribution of observed luminosity $\log_{10}(\nu
F_{\nu}D^{2})$ in the Spitzer IRAC and MIPS $24\mu$m bands vs MIR SED slope
$\alpha$ and concavity $a$ for all the YSOs in Gould’s Belt (first and third
columns) and the comparison with the similar distributions of the 200710 disk
models of Robitaille et al. (2006) (second and fourth columns; averaged over
the 50 model apertures), where $D$ is the average GAIA distance to each
observed protocluster and $D=250$ pc for the models. The red curve and error
bars show the average observed quantities (luminosity or concavity) and its
uncertainty in each $\alpha$ or luminosity bin, while the cyan error bars show
the sample dispersion in each bin. The meanings are similar for the model
curves (green) and error bars (green and orange; the green error bars are
usually very small) in the second and fourth column. The data of the curves
and error bars are given on the GitHub mentioned in Abstract.

## References

  * Adams et al. (1987) Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788, doi: 10.1086/164924
  * Alexander et al. (2006) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229, doi: 10.1111/j.1365-2966.2006.10294.x
  * Alibert et al. (2005) Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343, doi: 10.1051/0004-6361:20042032
  * Alves de Oliveira et al. (2012) Alves de Oliveira, C., Moraux, E., Bouvier, J., & Bouy, H. 2012, A&A, 539, A151, doi: 10.1051/0004-6361/201118230
  * André et al. (1993) André, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122, doi: 10.1086/172425
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
  * Azimlu et al. (2015) Azimlu, M., Martínez-Galarza, J. R., & Muench, A. A. 2015, AJ, 150, 95, doi: 10.1088/0004-6256/150/3/95
  * Bai (2016) Bai, X.-N. 2016, ApJ, 821, 80, doi: 10.3847/0004-637X/821/2/80
  * Baruteau et al. (2014) Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 667–689, doi: 10.2458/azu_uapress_9780816531240-ch029
  * Broekhoven-Fiene et al. (2014) Broekhoven-Fiene, H., Matthews, B. C., Harvey, P. M., et al. 2014, ApJ, 786, 37, doi: 10.1088/0004-637X/786/1/37
  * Burn et al. (2021) Burn, R., Schlecker, M., Mordasini, C., et al. 2021, A&A, 656, A72, doi: 10.1051/0004-6361/202140390
  * Caratti o Garatti et al. (2017) Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nature Physics, 13, 276, doi: 10.1038/nphys3942
  * Cieza et al. (2007) Cieza, L., Padgett, D. L., Stapelfeldt, K. R., et al. 2007, ApJ, 667, 308, doi: 10.1086/520698
  * Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485, doi: 10.1046/j.1365-8711.2001.04891.x
  * Contreras Peña et al. (2017) Contreras Peña, C., Lucas, P. W., Minniti, D., et al. 2017, MNRAS, 465, 3011, doi: 10.1093/mnras/stw2801
  * Crapsi et al. (2008) Crapsi, A., van Dishoeck, E. F., Hogerheijde, M. R., Pontoppidan, K. M., & Dullemond, C. P. 2008, A&A, 486, 245, doi: 10.1051/0004-6361:20078589
  * D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321, doi: 10.1086/320655
  * D’Alessio et al. (2006) D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314, doi: 10.1086/498861
  * Dodson-Robinson & Salyk (2011) Dodson-Robinson, S. E., & Salyk, C. 2011, ApJ, 738, 131, doi: 10.1088/0004-637X/738/2/131
  * Dong et al. (2015a) Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015a, ApJ, 809, L5, doi: 10.1088/2041-8205/809/1/L5
  * Dong et al. (2015b) Dong, R., Zhu, Z., & Whitney, B. 2015b, ApJ, 809, 93, doi: 10.1088/0004-637X/809/1/93
  * Dullemond & Dominik (2004a) Dullemond, C. P., & Dominik, C. 2004a, A&A, 421, 1075, doi: 10.1051/0004-6361:20040284
  * Dullemond & Dominik (2004b) —. 2004b, A&A, 417, 159, doi: 10.1051/0004-6361:20031768
  * Dullemond et al. (2001) Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957, doi: 10.1086/323057
  * Dullemond et al. (2007) Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 555, doi: 10.48550/arXiv.astro-ph/0602619
  * Dunham et al. (2015) Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11, doi: 10.1088/0067-0049/220/1/11
  * Emsenhuber et al. (2023) Emsenhuber, A., Burn, R., Weder, J., et al. 2023, A&A, 673, A78, doi: 10.1051/0004-6361/202244767
  * Emsenhuber et al. (2021) Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69, doi: 10.1051/0004-6361/202038553
  * Evans et al. (2009) Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321, doi: 10.1088/0067-0049/181/2/321
  * Fang et al. (2013) Fang, M., Kim, J. S., van Boekel, R., et al. 2013, ApJS, 207, 5, doi: 10.1088/0067-0049/207/1/5
  * Fang et al. (2009) Fang, M., van Boekel, R., Wang, W., et al. 2009, A&A, 504, 461, doi: 10.1051/0004-6361/200912468
  * Fischer et al. (2023) Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 355, doi: 10.48550/arXiv.2203.11257
  * Freedman & Diaconis (1981) Freedman, D., & Diaconis, P. 1981, Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 57, 453
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5, doi: 10.3847/0067-0049/224/1/5
  * Getman et al. (2014) Getman, K. V., Feigelson, E. D., & Kuhn, M. A. 2014, ApJ, 787, 109, doi: 10.1088/0004-637X/787/2/109
  * Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539, doi: 10.1088/0004-637X/690/2/1539
  * Greene et al. (1994) Greene, T. P., Wilking, B. A., Andre, P., Young, E. T., & Lada, C. J. 1994, ApJ, 434, 614, doi: 10.1086/174763
  * Guo et al. (2020) Guo, Z., Lucas, P. W., Contreras Peña, C., et al. 2020, MNRAS, 492, 294, doi: 10.1093/mnras/stz3374
  * Guo et al. (2024a) Guo, Z., Lucas, P. W., Kurtev, R., et al. 2024a, MNRAS, 528, 1769, doi: 10.1093/mnras/stad3700
  * Guo et al. (2024b) Guo, Z., Lucas, P. W., Kurtev, R. G., et al. 2024b, MNRAS, 529, L115, doi: 10.1093/mnrasl/slad201
  * Haisch et al. (2001) Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153, doi: 10.1086/320685
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Hartmann & Kenyon (1996) Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207, doi: 10.1146/annurev.astro.34.1.207
  * Harvey et al. (2007) Harvey, P., Merín, B., Huard, T. L., et al. 2007, ApJ, 663, 1149, doi: 10.1086/518646
  * Hsieh & Lai (2013) Hsieh, T.-H., & Lai, S.-P. 2013, ApJS, 205, 5, doi: 10.1088/0067-0049/205/1/5
  * Isella & Natta (2005) Isella, A., & Natta, A. 2005, A&A, 438, 899, doi: 10.1051/0004-6361:20052773
  * Jørgensen et al. (2009) Jørgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2009, A&A, 507, 861, doi: 10.1051/0004-6361/200912325
  * Karam & Sills (2023) Karam, J., & Sills, A. 2023, MNRAS, 521, 5557, doi: 10.1093/mnras/stad850
  * Kenyon et al. (1993) Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676, doi: 10.1086/173114
  * Kenyon et al. (1990) Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869, doi: 10.1086/115380
  * Kimura et al. (2016) Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257, doi: 10.1093/mnras/stw1531
  * Kirk et al. (2009) Kirk, J. M., Ward-Thompson, D., Di Francesco, J., et al. 2009, ApJS, 185, 198, doi: 10.1088/0067-0049/185/1/198
  * Kraus et al. (2012) Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19, doi: 10.1088/0004-637X/745/1/19
  * Krumholz (2011) Krumholz, M. R. 2011, in American Institute of Physics Conference Series, Vol. 1386, XV Special Courses at the National Observatory of Rio de Janeiro, ed. E. Telles, R. Dupke, & D. Lazzaro, 9–57, doi: 10.1063/1.3636038
  * Kuffmeier et al. (2018) Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642, doi: 10.1093/mnras/sty024
  * Kunitomo et al. (2021) Kunitomo, M., Ida, S., Takeuchi, T., et al. 2021, ApJ, 909, 109, doi: 10.3847/1538-4357/abdb2a
  * Kunitomo et al. (2020) Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-i. 2020, MNRAS, 492, 3849, doi: 10.1093/mnras/staa087
  * Lada (1987) Lada, C. J. 1987, in Star Forming Regions, ed. M. Peimbert & J. Jugaku, Vol. 115, 1
  * Larson (2003) Larson, R. B. 2003, Reports on Progress in Physics, 66, 1651, doi: 10.1088/0034-4885/66/10/R03
  * Liu et al. (2023) Liu, M., He, J., Ge, J., et al. 2023, ApJ, 943, 39, doi: 10.3847/1538-4357/aca658
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111, doi: 10.1088/0067-0049/186/1/111
  * Luhman et al. (2008) Luhman, K. L., Allen, L. E., Allen, P. R., et al. 2008, ApJ, 675, 1375, doi: 10.1086/527347
  * Mamajek (2009) Mamajek, E. E. 2009, in American Institute of Physics Conference Series, Vol. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii, 3–10, doi: 10.1063/1.3215910
  * Manara et al. (2023) Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 539, doi: 10.48550/arXiv.2203.09930
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192, doi: 10.1088/0004-6256/144/6/192
  * Mordasini et al. (2009) Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139, doi: 10.1051/0004-6361/200810301
  * Myers et al. (1987) Myers, P. C., Fuller, G. A., Mathieu, R. D., et al. 1987, ApJ, 319, 340, doi: 10.1086/165458
  * Palla & Stahler (2000) Palla, F., & Stahler, S. W. 2000, ApJ, 540, 255, doi: 10.1086/309312
  * Park et al. (2021) Park, W., Lee, J.-E., Contreras Peña, C., et al. 2021, ApJ, 920, 132, doi: 10.3847/1538-4357/ac1745
  * Pignatale et al. (2018) Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E. 2018, ApJ, 867, L23, doi: 10.3847/2041-8213/aaeb22
  * Pokhrel et al. (2020) Pokhrel, R., Gutermuth, R. A., Betti, S. K., et al. 2020, ApJ, 896, 60, doi: 10.3847/1538-4357/ab92a2
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Raymond & Morbidelli (2022) Raymond, S. N., & Morbidelli, A. 2022, in Astrophysics and Space Science Library, Vol. 466, Demographics of Exoplanetary Systems, Lecture Notes of the 3rd Advanced School on Exoplanetary Science, ed. K. Biazzo, V. Bozza, L. Mancini, & A. Sozzetti, 3–82, doi: 10.1007/978-3-030-88124-5_1
  * Rebull et al. (2010) Rebull, L. M., Padgett, D. L., McCabe, C. E., et al. 2010, ApJS, 186, 259, doi: 10.1088/0067-0049/186/2/259
  * Robitaille (2008) Robitaille, T. P. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 387, Massive Star Formation: Observations Confront Theory, ed. H. Beuther, H. Linz, & T. Henning, 290, doi: 10.48550/arXiv.0711.4369
  * Robitaille et al. (2006) Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256, doi: 10.1086/508424
  * Rosotti et al. (2013) Rosotti, G. P., Ercolano, B., Owen, J. E., & Armitage, P. J. 2013, MNRAS, 430, 1392, doi: 10.1093/mnras/sts725
  * Schib et al. (2021) Schib, O., Mordasini, C., Wenger, N., Marleau, G. D., & Helled, R. 2021, A&A, 645, A43, doi: 10.1051/0004-6361/202039154
  * Shu et al. (1987) Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23, doi: 10.1146/annurev.aa.25.090187.000323
  * Soderblom et al. (2014) Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 219–241, doi: 10.2458/azu_uapress_9780816531240-ch010
  * Varnière et al. (2006) Varnière, P., Blackman, E. G., Frank, A., & Quillen, A. C. 2006, ApJ, 640, 1110, doi: 10.1086/498933
  * Vorobyov & Basu (2015) Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115, doi: 10.1088/0004-637X/805/2/115
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1, doi: 10.1086/422992
  * Whitney & Hartmann (1993) Whitney, B. A., & Hartmann, L. 1993, ApJ, 402, 605, doi: 10.1086/172163
  * Whitney et al. (2003) Whitney, B. A., Wood, K., Bjorkman, J. E., & Cohen, M. 2003, ApJ, 598, 1079, doi: 10.1086/379068
  * Wilking et al. (1989) Wilking, B. A., Lada, C. J., & Young, E. T. 1989, ApJ, 340, 823, doi: 10.1086/167439
  * Wood et al. (2002) Wood, K., Lada, C. J., Bjorkman, J. E., et al. 2002, ApJ, 567, 1183, doi: 10.1086/338662
  * Zhu et al. (2010) Zhu, Z., Hartmann, L., & Gammie, C. 2010, ApJ, 713, 1143, doi: 10.1088/0004-637X/713/2/1143