© The Authors 2026

1 Introduction

Constraining circumstellar disk lifetimes is crucial to understanding the timescale for planet formation (Morbidelli & Raymond 2016). Circumstellar disks provide the material from which planets form, but they eventually disperse due to accretion onto and radiation from the central star. External influences, such as nearby massive stars, can further impact disk dispersal and lifetime. These disks originate during the star formation process, as contracting cores in molecular clouds establish a rotating disk around the forming star, supplying additional material.

Haisch et al. (2001) pioneered disk lifetime estimation by analyzing the fraction of young stars with infrared (IR) excess, indicative of circumstellar disks, across stellar clusters of varying ages (see also the review by Hillenbrand 2005). Their linear fit yielded a disk lifetime of about 6 Myr, defined by the x-axis intercept in the disk fraction versus stellar cluster age plot. Hernández et al. (2007) collected for the first time the disk fractions and ages for multiple young stellar clusters in a systematic way. Subsequent studies to the first timescale estimate, employing more realistic exponential decay models, broadened the estimated range. Typical values of the e-folding time now fall between 2 and 4 Myr (Fedele et al. 2010; Ribas et al. 2014, 2015; Richert et al. 2018) and around 8 Myr (Michel et al. 2021). Although some studies use alternative definitions, the results remain broadly consistent. For instance, Hernández et al. (2008) identified a rapid disk fraction decay at approximately 5 Myr, tracing the complete dispersion of disks rather than e-folding timescales. Pfalzner et al. (2022) found a median lifetime of 5–10 Myr using exponential fits on the relation of disk fraction versus stellar cluster age and Gaussian fits to the lifetime distribution of disks (for an explanation of this, see Pfalzner 2022). It is important to note that all these estimates are susceptible to selection effects (Pfalzner et al. 2014, 2022) and stellar mass dependencies (Carpenter et al. 2006; Roccatagliata et al. 2011; Luhman & Mamajek 2012; Yasui et al. 2014; Ribas et al. 2015).

The determination of disk lifetimes is subject to signifi-cant uncertainty due to several observational and theoretical challenges. One of the primary sources of uncertainty lies in the estimation of stellar ages, which are typically derived from Hertzsprung-Russell (HR) diagram positions and theoretical pre-main sequence (PMS) evolutionary models. The situation is improved when considering the age distribution of a population of stars, but even in this case, models exhibit systematic discrepancies, as shown, for example, in Ratzenböck et al. (2023b).

Using populations of coeval stars to mitigate age estimation uncertainties introduces its own set of challenges. One significant issue is that these populations are at different distances from Earth, which usually implies different completeness across a sample. More distant populations are harder to observe, leading to an incomplete census of stars, especially for the fainter, but numerous, low-mass end. This distance-related bias can skew the age distribution and disk lifetime estimates, as the observed sample might not be representative of the true population. To make matters all the more challenging, different young stellar populations might have been exposed to different physical conditions, such as the presence of disk-destructive high-UV environments around massive stars, introducing additional uncertainties in a comparison between them.

Finally, the limited availability of well-characterized samples of young stars in the 5–20 Myr age range exacerbates these uncertainties, making it challenging to derive a robust statistical distribution of disk lifetimes. Consequently, while broad trends in disk dissipation are evident, precise quantification of disk lifetimes and the dominant mechanisms driving their evolution remain areas of active research. For instance, Richert et al. (2018) conducted one of the largest studies on the lifetimes of inner dust disks using X-ray and IR photometry for 69 young stellar clusters in 32 nearby star-forming regions, only covering ages less than or equal to 5 Myr.

In this study, we take an approach designed to reduce key uncertainties in disk lifetime estimates (age and completeness) by focusing on a single star-forming region. Unlike previous studies that combined multiple young stellar clusters, our analysis is based on a homogeneous sample within a single star-forming region, the nearby Scorpius-Centaurus OB association (Sco-Cen). Ratzenböck et al. (2023a) introduced a novel clustering method, significance mode analysis (SigMA), which they first applied to Gaia DR3 data (Gaia Collaboration 2023) for Sco-Cen. The authors identified 34 co-spatial, co-moving, and coeval stellar clusters¹ within Sco-Cen, as also confirmed in Ratzenböck et al. (2023b). By deriving isochronal ages from well-defined isochrones in HR diagram, they validated their cluster selection and found ages ranging from 3 to 21 Myr. These ages are based on PARSEC model isochrones (see the discussion on ages from different models in Ratzenböck et al. 2023b). This clustering approach mitigates the issue of age mixing in a single star-forming region, enabling a coherent study of cluster properties over time within a single stellar association at a single distance, minimizing differences in the completeness of the various clusters.

The goal of this work is to construct the first disk fraction versus stellar cluster age plot for a single, nearby, and well-defined stellar population. With this homogeneous sample covering a broad age range, we can refine the interpretation of disk evolution as a function of stellar cluster age. A key advantage of this approach is the consistent determination of both cluster membership and age, ensuring a uniform analysis across all clusters.

2 Data

Our analysis is based on the set of stellar clusters identified within the Sco-Cen OB association by Ratzenböck et al. (2023a), who reported 37 candidate clusters. These stellar clusters have been selected as co-spatial and co-moving statistical overdensi-ties in the five-dimensional phase space (three positional coordinates and two velocity components). Subsequently, Ratzenböck et al. (2023b) found that 34 of these are more likely to be physically associated with Sco-Cen. From these, we excluded the cluster Centaurus-Far from our sample due to its potentially high level of contamination from older field stars, which renders its age estimate unreliable. The remaining 33 clusters, comprising a total of 12 873 stellar sources, form the basis of our study. Ratzenböck et al. (2023b) derived isochronal ages and their uncertainties using Gaia photometry and PARSEC stellar evolutionary models (Marigo et al. 2017). Given this approach, we are not including the embedded population in our sample and these are also not accounted for in the age estimation. The age uncertainties are only used for visualization in our work and are not included in the inference process described in Sect. 3. The cluster ages span a range of approximately 3–21 Myr, covering the critical period during which most circumstellar disks are expected to dissipate. For each source, we use the Gaia source_id, cluster membership label, and corresponding age. For each stellar source, we adopted the age and age uncertainty of its parent cluster.

To assess the presence of circumstellar disks, we identified IR excess via cross-matching with near- and mid-infrared (NIR and MIR) photometric catalogs. Infrared measurements were obtained from the AllWISE catalog (Cutri et al. 2013), using the Gaia–AllWISE best-neighbor cross-match. We found valid WISE counterparts for 10 087 sources, corresponding to 78.4% of our initial sample.

In addition to MIR data, we incorporated NIR measurements from the Two Micron All Sky Survey (2MASS) (Skrutskie et al. 2006). Since 2MASS data are already integrated into the All-WISE catalog for sources with MIR detections, this subset is sufficient for our analysis. All disk selection methods described below rely on a combination of WISE and 2MASS photometry to identify IR excess indicative of circumstellar material.

3 Methods

3.1 Detecting circumstellar disks

We identified circumstellar disks by detecting IR excess in young stellar objects (YSOs) across all stellar clusters analyzed in Sco-Cen. Throughout this work, we refer to sources exhibiting IR excess as disk-bearing (DB) and those without as disk-less (DL). While no further subclassification into YSO classes is applied, the MIR colors of most disk-bearing sources in our Gaia-selected sample suggest they are predominantly Class II candidates.

Our primary method for disk identification relies on a combination of IR color–color diagrams (CCDs), which serve as our main disk selection approach. As a complementary method, we also applied a selection based solely on the shape of the (extinction-uncorrected) IR spectral energy distribution (SED). In both cases, we aim to identify sources with clear IR excess indicative of the presence of a circumstellar disk.

These two approaches, referred to as the CCD selection and SED selection, respectively, were applied independently to compute disk fractions for each stellar cluster. By comparing the resulting disk lifetimes from each method, we were able to estimate the level of systematic uncertainty introduced by the choice of disk selection strategy. A detailed description and comparison of both selection methods are provided in Appendix A.

In addition, we compared our selections to that of Luhman (2022), who analyzed disk-bearing sources in a subregion of Sco-Cen encompassing 25 of the 33 stellar clusters included in our study. This externally derived selection of disks is referred to as the Luhman selection.

Figure 1 shows the sky distribution of the disk-bearing sources (identified via the CCD selection) overlaid on the diskless population, as discussed further in Sect. 4.

Fig. 1 Plot of galactic longitude versus galactic latitude showing the region of Sco-Cen. The sources are shown based on their selection, blue circles for disk-bearing, and orange crosses for disk-less sources.

3.2 Fitting disk fraction versus age

We aim to derive the disk decay parameter from the empirically found functional relationship between a stellar cluster’s age and its average disk fraction. With sources separated into disk-bearing and disk-less, the disk fraction, f_i, for each stellar cluster, i, is calculated as $$f_i=\frac{\#D{B}_i}{\#D{B}_i+\#D{L}_i},$$(1) where #DB_i is the number of disk-bearing sources in the stellar cluster i and #DL_i is the number of disk-less sources. We choose the function here on the basis of two factors. First, the observed shape of the relation between the disk fraction and the age t (see Fig. 2) empirically resembles an exponential decay. Second, exponential decay was already suggested to be an appropriate solution in the earlier literature (Fedele et al. 2010; Ribas et al. 2014, 2015; Richert et al. 2018; Michel et al. 2021; Pfalzner et al. 2022). We modeled the disk fraction as a function time using the following relationship, $$f(t)=f_0\cdot e^{-\frac{t-t_0}{\tau }}.$$(2)

Here, t (Myr) denotes the age of a stellar cluster, and τ (Myr) is the characteristic decay parameter, representing the decay time of the disks and the main parameter we are interested in. Subtracting t₀ (Myr), referred to as the shift parameter, in the exponential’s numerator shifts the function horizontally along the age axis t (the x-axis). Introducing t₀ allows for different onset times of the exponential decay to be modeled directly, accounting for situations in which disks do not dissipate immediately at the birth of stellar clusters. The parameter f₀ denotes the intercept or initial disk fraction at the start of the exponential decay time t₀. As the function is shifted in a parallel manner along the x-axis, this point is evaluated at t − t₀ = 0.

3.3 Parameter inference

This work aims to constrain the decay parameter τ using disk fractions and ages of stellar clusters in Sco-Cen. To this end, we constructed a likelihood function based on the parametric disk evolution model in Eq. (2) and used MCMC sampling to estimate the posterior distribution of τ. This approach offers a principled way to propagate observational uncertainties and quantify credible intervals on τ and other parameters of interest. We provide a detailed discussion of our likelihood function in Appendix B. In a regime with limited data, the choice of priors becomes particularly important, as it can notably influence the inferred posterior distributions. Whenever possible, we adopted weakly informative priors to regularize the inference and prevent overfitting, while still allowing the data to drive the results. These priors, representing our “fiducial model”, were selected to reflect plausible values for astrophysical parameters based on previous studies, without imposing overly strong constraints. In the following, we briefly highlight these choices and refer to Appendix B for a more in-depth discussion.

3.3.1 Fiducial model

We used the prior probability distributions that encode our physical expectations while allowing the data to drive the final results. For the decay timescale, τ, we aimed to ensure positive values while remaining flexible across a wide range of physically plausible timescales. The half-Cauchy distribution is a good candidate to encode these prior considerations, as its heavy tail imposes minimal constraints on large values, allowing the data to determine whether disk dispersal is rapid or gradual. For the initial disk fraction or intercept f₀, we use a Beta distribution centered at 0.83, reflecting the findings by Pfalzner & Dincer (2024) that initial disk fractions likely fall between 0.65 to 1.0. The Beta distribution is ideal for modeling fractions because it is naturally bounded between 0 and 1 and its shape can be adjusted to reflect varying degrees of prior confidence. For the shift parameter, t₀, we adopted a half-Normal distribution, which assigns a decreasing probability to large positive values. This choice reflects our physical expectation that disk dispersal likely begins close to the cluster formation epoch rather than being substantially delayed, while still permitting larger offsets if strongly indicated by the data. The half-Normal’s relatively rapid decay towards zero provides more structure than the heavy-tailed half-Cauchy, which is appropriate for t₀ since we have stronger prior physical intuition about this parameter than about the decay timescale, τ. Complete specifications of all prior distributions are provided in Appendix B.

3.3.2 Robustness checks and model variants

To assess the robustness of our results and evaluate the sensitivity of the decay time, τ, to prior choices and parameter constraints, we also considered two alternative models with distinct assumptions. First, we implemented a “data-driven” model with fully uninformative priors (in contrast to the physically constrained priors of the fiducial model), allowing for a broad exploration of the parameter space, including potentially unphysical values (e.g., an intercept greater than 1). This model yields a posterior distribution that is effectively driven by the likelihood function alone, enabling an assessment of the information content in the data independent of prior regularization.

Second, we introduced a model that omits the shift parameter, t₀, referred to as the “no-shift” model. This approach, commonly used in the literature (Fedele et al. 2010; Richert et al. 2018; Michel et al. 2021; Pfalzner et al. 2022), provides a direct comparison to previous studies and serves as a test to determine whether incorporating t₀ significantly affects the inferred disk decay time.

To ensure that our conclusions are not biased by the choice of disk selection method, we applied the fiducial model across the CCD selection, the SED selection, and the Luhman selection of disks from Luhman (2022). Comparing results across these selections allows us to assess the stability of the inferred decay timescale, τ, and identify potential systematics associated with disk selection criteria. A detailed description of each selection method and its implementation is provided in Appendix A.

Additionally, we applied the fiducial model on two further subsamples by splitting the CCD selection into a possible binary and single star sample to study the influence of binary stars. Kraus et al. (2012) showed that binarity could decrease the measured disk fraction. To identify binary candidates, we used the Gaia RUWE parameter (Penoyre et al. 2022a,b; Castro-Ginard et al. 2024), while a rough quality selection of RUWE ≥ 1.4 indicates Gaia observed sources, which are more likely multiple stellar systems. Castro-Ginard et al. (2024) showed that the threshold on the RUWE parameter varies with the position in the sky. However, we adopted a threshold of 1.4, consistent with Ratzenböck et al. (2023b), who also used this quality criterion to exclude binary candidates in the Sco-Cen sample.

Fig. 2 Disk fraction versus age. The black data points represent the stellar clusters, on the left for the CCD selection, in the middle for the SED selection and on the right for the Luhman selection. The Luhman selection only includes 25 out of 33 stellar clusters because the study by Luhman (2022) encompasses a smaller region compared to Ratzenböck et al. (2023a). Shown in green-blue in the right plot are the missing stellar clusters (from the CCD selection). In the left plot, the red line represents the median of the fiducial+CCD result, which is also shown in the middle and right plot with reduced transparency. The orange dash-dotted line in the middle plot gives the median of the fiducial+SED result. The blue dotted line is the median of the fiducial+Luhman result. The gray area shows the 1σ confidence interval for each median result. The uncertainties in ages are from Ratzenböck et al. (2023b) and the disk fraction uncertainties are derived from the underlying counting process which can be modeled as a simple Bernoulli trial.

4 Results

4.1 Disk fraction per stellar cluster

We selected 10 016 sources showing the presence or absence of a circumstellar disk in Sco-Cen, resulting from the analysis of sources characterized by IR excess using the CCD selection, as outlined in Appendix A. Using the SED selection, we assessed the presence or absence of a disk for 10 021 sources, with 10 008 appearing in both selections. The Luhman selection contains in total 8137 sources. We used the CCD selection as our default disk selection in this work, since the various band combinations allow a more thorough identification of disk-bearing sources and it is a comparable method to some approaches from the literature (Koenig & Leisawitz 2014; Großschedl et al. 2019, 2021; Luhman 2022). The CCD selection excludes stellar sources that are reddened due to extinction, which is why we omitted the extinction correction. We refer to Appendix A.1 for more details.

Figure 1 shows the sky distribution of the disk-bearing young stars from the CCD selection on top of the remaining diskless sources. This figure highlights the extent of Sco-Cen as selected in Ratzenböck et al. (2023a). Disk-bearing sources are distributed throughout the entire Sco-Cen region, with several noticeable over-densities. In particular, disks are more concentrated on top of known star-forming regions, which are still containing molecular clouds, mostly at the periphery of the association.

The resulting disk statistics per stellar cluster are given in the Table C.1, including the names and ages of each stellar cluster. We list the number of disk-bearing and disk-less sources, and the disk fraction for each of our selection approaches. These disk fractions are used in the following section to determine the disk decay time.

4.2 Derived decay times

Figure 2 presents the relationship between disk fraction and stellar cluster age for all stellar clusters in Sco-Cen, with results shown separately for the three disk selection methods. The left, middle, and right panels correspond to disk fractions derived using the CCD, SED, and Luhman selections, respectively. Since the underlying cluster ages remain fixed across all panels, differences between the subplots arise solely from variations in disk selection methods, namely, in the y-axis values. The stellar clusters not included in the Luhman selection are shown in green-blue (from the CCD selection).

In each panel, we overlay the median posterior prediction from the fiducial model fitted to the corresponding disk selection, along with the 1σ credible interval indicated by the shaded gray region. The exponential function is evaluated for each posterior sample of the parameter triplet – decay time τ, intercept f₀, and shift t₀ – starting from 0 Myr. For ages t < t₀, the model predicts a constant disk fraction equal to the intercept, resulting in an initial plateau followed by an exponential decline once t > t₀. The median of the model ensemble yields a smooth transition between these regimes. Age uncertainties are adopted from Ratzenböck et al. (2023b), while uncertainties in the disk fraction arise from the underlying source selection process. Since each star is either disk-bearing or disk-less, the disk count in each cluster follows a Bernoulli process, and the corresponding uncertainty can be modeled using a Beta-Binomial distribution (see Appendix B). Overall, we found that the exponential decay model provides a good empirical description of the disk fraction across the observed age range.

Figure 3 compares the inferred decay times τ across the different model variants and disk selection methods, as sum-marized in Table 1. To facilitate comparison, we display the posterior probability distributions as violin plots, which illustrate both the distribution shape and credible intervals. Overall, the decay times obtained from the various models are consistent within their respective uncertainties. However, the alternative disk selection methods (SED and Luhman) tend to yield slightly larger median decay times. This suggests that the choice of the disk selection approach (which can systematically bias disk fraction estimates) has a measurable impact on the inferred decay timescale.

Table 1 summarizes the detailed results of our exponential fits to the disk fraction as a function of stellar cluster age. For each combination of the model variant (fiducial, data-driven, no-shift) and disk selection method (CCD, SED, Luhman), we report the median values of the fitted parameters (τ, f₀, t₀), along with their 1σ uncertainties, defined via the 68% highest density interval (HDI). The CCD selection serves as our primary reference sample.

In the final two rows of the table, we present results from two subsamples of the CCD selection, separated on the basis of stellar multiplicity using the Gaia RUWE parameter. Sources with RUWE < 1.4 are more likely to be single stars, while those with RUWE ≥ 1.4 are candidate multiple systems. For the single-star subsample, the posterior distribution of the decay time, τ, shifts toward higher values, placing it between the results of the SED and Luhman selections. In contrast, the multiple system candidate subsample yields a lower τ value. Notably, the decay time inferred from the full CCD sample using the fidu-cial model lies between those of the two subsamples, suggesting that stellar multiplicity might influence disk lifetimes, by slightly underestimating the lifetime.

Among the tested model and selection combinations, we identified the fiducial model paired with the CCD-based disk selection as the most robust and reliable configuration. Compared to the SED and Luhman approaches, the CCD method incorporates a broader set of photometric bands and color combinations, providing a more comprehensive basis for disk selection methods. The fiducial model also encodes the most physical prior information, although its influence on the posterior estimates is modest. Based on this preferred configuration, we adopted a characteristic disk decay time of 5.8 ± 0.3 Myr as our main result.

Taken together, the full range of tested models and data suggests that the characteristic disk dissipation timescale in Sco-Cen likely falls between approximately 5 and 7 Myr. To our knowledge, this represents the first such determination within a single, coherent star-forming region. We emphasize that our results are specific to the Sco-Cen association and should not be assumed to generalize across star-forming regions without further investigation.

Fig. 3 Violin plot showing the posterior probability density functions (PDFs) of the disk decay time parameter τ. The x-axis denotes the decay time in Myr, while the y-axis labels each model and disk selection method. The number N displayed next to each distribution indicates the number of sources used in the corresponding fit; this number remains constant for the top three rows, where only the model configuration varies. Posterior distributions are grouped and color-coded: red indicates variations in prior assumptions (with the fiducial model highlighted at the top), green represents the supplementary disk selection methods (SED, Luhman) and blue corresponds to subsamples based on stellar multiplicity (single vs. candidate binaries). See Table 1 for the corresponding numerical results.

5 Discussion

The main result of this work is the precise derivation of disk lifetime for a single star-forming region. We found a local disk lifetime in the Sco-Cen OB association to be 5.8 ± 0.3 Myr, approximately a factor of two longer than most previous estimates, suggesting that planet formation might take more time than previously assumed. The most straightforward explanation for our longer lifetime estimate is that in this work we were able to sample, for the first time, the (essentially) complete stellar mass spectrum down to the hydrogen-burning limit for a single stellar population of young stars located at the same distance, minimizing systematic uncertainties. As previously argued in Pfalzner et al. (2022), earlier estimates of short disk lifetimes might be biased. This potential bias arises from combining disk samples that include clusters at different distances, leading to an overrepresentation of high-mass stars. Since high-mass stars tend to disperse their disks at earlier stages, the resulting disk lifetimes are not representative of the broader stellar population, which is dominated by low-mass stars. However, we calculated the disk decay time for the entire mass sample in our data, which can end up having an influence on our result (Carpenter et al. 2006; Roccatagliata et al. 2011; Luhman & Mamajek 2012; Yasui et al. 2014; Ribas et al. 2015).

5.1 Model and selection effects on disk decay times

The goal of this work is to constrain the characteristic decay time of circumstellar disks. Figure 3 compares the posterior distributions of the decay time parameter τ across different model configurations and disk selection methods. Overall, the decay times fall within a range of approximately 5–7 Myr. The fidu-cial+CCD result shows a mostly symmetric posterior centered at ∼5.8 Myr, while the data-driven+CCD model yields a broader and more asymmetric distribution, with a longer tail toward higher values. Despite these differences in shape, all three models tested with the CCD selection (fiducial, data-driven, and no-shift) produce consistent median decay times (see Table 1), suggesting the results are robust against modeling choices.

We examined the impact of the disk selection method on the inferred decay time. Compared to the CCD selection, both the SED- and Luhman-based selections yield longer decay times and lower intercepts. These differences arise from the way each method selects disks and from differences in the underlying stellar samples. In particular, the fiducial+SED result appears to underrepresent younger clusters with high disk fractions, resulting in a flatter fit and thus a longer inferred decay time. We did not correct the SEDs for extinction as we considered it a second-order effect when influencing the selection of disk-bearing or disk-less sources; rather, the cut at one specific value of −2 can have an influence on the disk fractions. Similarly, the Luhman selection omits several young clusters (e.g., B59, Chamaeleon-1/2), leading to a lower intercept and an even flatter decay curve. The missing clusters are predominantly younger and have higher disk fractions, which steepens the CCD-based decay curve relative to the Luhman-based result. These trends are clearly visible in Fig. 2 and quantitatively reflected in Table 1.

To explore another possible source of systematic variation, we also divided the CCD sample into candidate single and multiple systems using the Gaia RUWE parameter. While the resulting decay times differ (consistent with earlier suggestions that binarity reduces disk lifetimes), the interpretation remains tentative. The shorter decay time for the high-RUWE subsample might reflect either a genuine multiplicity effect or the influence of small-number statistics, given the limited size of the binary candidate sample (see the number N in Fig. 3).

Despite these differences, all selection methods exhibit a clear exponential decline in disk fraction with age. The uncertainty in the decay timescale is lowest at older ages, where the data are densest and the fits converge. At younger ages, differences in disk selection and sample completeness become more significant, leading to larger variation in the inferred decay time. Still, the shift parameter t₀ remains consistent across all selection methods, and the general trend supports a characteristic exponential disk dissipation timescale for Sco-Cen.

5.2 Comparison with Pfalzner & Dincer 2024

Traditionally, the Sco-Cen OB association has been separated into three substructures: Upper-Scorpius, Upper-Centaurus-Lupus, and Lower-Centaurus-Crux (Blaauw 1964). However, we now know from Gaia, and in particular from the work of Ratzenböck et al. (2023b) that it contains more than 30 stellar clusters. The earlier and coarser subdivisions contributed to mixing of subpopulations of different ages and different age estimations, especially when using data that predated Gaia (see discussion on the Upper-Sco age, pre- and post-Gaia in Ratzenböck et al. 2023b).

Next, we compared our results with the results in Pfalzner & Dincer (2024), one of the most recent work on the topic. These authors fit a disk lifetime distribution to a set of stellar clusters compiled from various literature sources in the local Milky Way. Their sample includes clusters and disk fractions obtained through different methods. They focus on stars with spectral types from M3.7 to M6 to limit the mass range and test several probability distributions. In their study, Sco-Cen appears as two sub-populations: Upper Scorpius and Upper-Centaurus-Lupus/Lower-Centaurus-Crux. Their results yield median disk lifetimes on the order of 5–10 Myr, comparable to our findings. Still, and for comparison, what is typically called Upper Scorpius in the literature corresponds to about ten clusters in our sample, with ages ranging from 3 up to 19 Myr, allowing for a more precise determination of disk lifetimes.

In Pfalzner & Dincer (2024), a cut in spectral type (or mass) further refines the sample and excludes the effect of mass-dependent disk lifetimes, which are believed to be shorter for higher-mass stars (Carpenter et al. 2006; Roccatagliata et al. 2011; Luhman & Mamajek 2012; Yasui et al. 2014; Ribas et al. 2015). We do not account for stellar mass, so our derived decay times might be reduced by the inclusion of sources spanning a wide range of masses. Pfalzner & Dincer (2024) observe consistently large standard deviations of around 6 Myr, suggesting very broad lifetime distributions. Our derived disk decay times, shown in Fig. 3, also cover a broad range, but the standard deviations in Table 1 are substantially smaller. This discrepancy could arise because we apply an exponential function and examine the distribution of its parameters.

Pfalzner & Dincer (2024) also explored different initial disk fraction values. Their best fit occurs near an initial disk fraction of 0.65, but we usually find larger intercepts, f₀, except when using data from Luhman (2022) or the RUWE < 1.4 model. It is important to note that our intercept is defined at t − t₀ = 0, so an exponential function starting at 0.65 would already lie below the disk fraction of some stellar clusters.

Additionally, Pfalzner & Dincer (2024) show that different functional forms might be needed to fit these properties. Our data, however, are well-represented by an exponential decay over the age range we cover, and both the shifted and unshifted exponential models work. Due to limited data below 3 Myr, we cannot constrain the shape of the disk fraction–age relationship in very young systems. Nonetheless, our inference pipeline suggests an intercept f₀ that is not equal to 1 (see also the discussion below).

Finally, our approach excludes younger stellar clusters that have been identified with different methods or that are located in different regions at different distances. This consistency in selection, achieved with SigMA, reduces systematic uncertainties. Clusters in different regions might have, for example, different formation conditions, different levels of UV flux, so merging them could lead to misleading results. Our approach thus benefits from a homogeneously selected and analyzed sample, both in the clustering process and in the disk selection methods.

5.3 The initial disk fraction, f₀

Richert et al. (2018) and Michel et al. (2021) explored the effect of different initial disk fractions and find that this value might be below 1 for multiple reasons, such as binarity or disks dissipating rapidly or not forming at all. In our analysis, we did not fix the numerical value of the intercept f₀ at t − t₀ = 0; rather, we left the intercept as a free parameter to be inferred. This allows for the data to inform the posterior distribution, while incorporating prior knowledge that f₀ must lie within the physical range [0, 1].

Our results consistently favor initial disk fractions below unity, with posterior values typically spanning 0.6–0.85. This is in line with findings by Michel et al. (2021), who report f₀ ≈ 0.65 and attribute the reduction to binarity. However, we find that even in the RUWE < 1.4 subsample – expected to contain fewer unresolved binaries – the inferred f₀ remains below that of the full fiducial+CCD sample. This suggests that binarity alone might not account for the sub-unity intercept and that other mechanisms, such as intrinsically low initial disk fractions or rapid early disk dispersal, could also play a role as proposed by Richert et al. (2018).

It is important to note that our data do not extend below 3 Myr and, thus, we cannot empirically constrain the behavior of the disk fraction at very young ages. As a result, the intercept, f₀, is influenced by a transition from a data-rich regime (older than ∼5 Myr) to a data-poor regime at younger ages. In this early-age regime, the posterior distribution is increasingly shaped by the prior. Therefore, while our results support a sub-unity initial disk fraction, further observations of younger stellar populations will be needed to confirm the early-time behavior of disk evolution.

5.4 The effect of stellar ages on inferred decay times

Our analysis has not explicitly propagated the uncertainties in individual cluster ages through the inference process. The ages employed here are maximum a posteriori (MAP) estimates reported in Ratzenböck et al. (2023b) derived from PARSEC isochrones using BP-RP. Isochronal age estimates are inherently subject to systematic uncertainties that depend on input physics and stellar evolution models. The application of identical models across the full range from very young (~1 Myr) to older (~20 Myr) clusters could introduce systematic biases (Bell et al. 2013). Furthermore, the neglect of accretion processes in pre-main sequence models can lead to overestimated ages for massive stars (Hosokawa et al. 2011). In the following, we apply our inference pipeline to multiple age estimates across different model families to estimate the systematic uncertainty of our reported e-folding time.

To assess the systematic impact of age uncertainty on our results, we re-analyze our fiducial model using alternative age estimates from Ratzenböck et al. (2023b) and consider other reported literature ages from Kerr et al. (2021), combined with the CCD selection. As demonstrated by Ratzenböck et al. (2023b), one source of systematic uncertainty stems from the choice of photometric color space. We therefore apply our methodology to age sets derived from two isochrone model families, PARSEC (Bressan et al. 2012) and Baraffe (Baraffe et al. 2015), using two different Gaia color combinations: BP-RP and G-RP.

When adopting PARSEC-G-RP ages, we derive a median decay time of ${6.09}_{-0.41}^{+0.32}$ Myr, which agrees with our primary result within uncertainties. The Baraffe-G-RP combination yields ${4.46}_{-0.27}^{+0.25}$ Myr, somewhat lower than our fiducial estimate. This discrepancy is significant and demonstrates the sensitivity to the age determination methods employed. Ratzenböck et al. (2023b) showed that Baraffe-G-RP ages are consistent, within uncertainties, with ages determined from PARSEC isochrone fits.

The Baraffe-BP-RP age combination exhibits the strongest systematic offset, consistently underestimating ages relative to the PARSEC-BP-RP reference. This translates to a decay time of ${3.28}_{-0.24}^{+0.30}$ Myr, which approaches the shorter timescales reported in the literature. Despite this apparent consistency with previous studies, Ratzenböck et al. (2023b) demonstrate that this age estimate represents a significant outlier compared to all other age approximations.

For additional perspective, we compare our results using ages from Kerr et al. (2021), which provide an independent cluster analysis of the Sco-Cen region. Ratzenböck et al. (2023a) provided a detailed comparison to Kerr et al. (2021) and we selected groups clearly identifiable in both analyses as an additional age estimate for comparison. Applying our methodology to the subset of clusters with cross-matched ages, determined in Kerr et al. (2021), yields a decay time of ${6.05}_{-0.56}^{+0.42}$ Myr, further demonstrating the sensitivity of our inferred parameters to the adopted age scale.

The systematic variations demonstrated here underscore that age determination represents a fundamental limitation in constraining disk evolution timescales and likely contributes significantly to the observed discrepancies between different studies in the literature. Taking this into consideration, we expand the range of decay time values to 4–7 Myr with a skewness towards larger values. Similar to Ratzenböck et al. (2023b) considering the Baraffe-BP-RP ages as outlier, we see the derived decay time using these ages as one. We show the disk fraction versus age plot with the inferred functional relation same as Fig. 2 in the Appendix in Fig. C.4. In a recent work, Fang & Herczeg (2025) use the disks from Luhman (2022) to fit a standard exponential decay to the disk fraction versus age relation using the stellar clusters from Ratzenböck et al. (2023a) and obtain a decay time of about 6 Myr for K- and M-type stars. They find a spread of decay times depending on the age determination method with even larger values when accounting for stellar cool spot coverage in the age determination.

5.5 Parameters and their influence

We briefly examine how the model parameters, decay time τ, intercept f₀, and shift t₀, interact and influence the shape of the fitted disk fraction curve. The shift parameter t₀ allows for a delay in the onset of exponential decay, which can introduce degeneracies with the intercept f₀. This is evident in the corner plots (Appendix C), where a strong correlation between f₀ and t₀ is visible, especially in the fiducial and data-driven models.

Excluding the shift parameter (in the no-shift model) leads to a stronger correlation between τ and f₀, as the model compensates for the fixed onset by adjusting the slope and initial value. Including t₀ helps to decouple these parameters and improves model flexibility, though t₀ itself remains weakly constrained in most cases. Despite these correlations, all three model variants (fiducial, data-driven, no-shift) yield consistent decay timescales, indicating that the main conclusions are robust to the inclusion or exclusion of the shift parameter.

5.6 Interpretation, limitations, and context

First, our sample is based on Gaia-selected stellar members, which excludes highly extincted or embedded sources. This limitation likely results in an underestimation of disk fractions, particularly in younger clusters still embedded in molecular clouds. If these obscured sources were included, the initial disk fractions would likely be higher. For a given disk fraction at older ages, this would require a steeper exponential decline, thus a shorter decay timescale, to connect the higher starting point to the same endpoint. Therefore, our exclusion of embedded sources could lead to a slight overestimation of the disk lifetime. Nonetheless, extinction toward Sco-Cen is relatively low overall, which helps mitigate this bias and makes the region particularly well-suited for this type of analysis.

Second, we note that the wavelength range used for disk selection can influence the measured disk lifetime. Ribas et al. (2014) find that the MIR data yield systematically longer disk lifetimes (compared to the NIR) due to improved sensitivity to IR excess. Our use of WISE MIR bands (W3, W4) in both the CCD and SED-based selection methods helps mitigate this concern by incorporating wavelength regimes more sensitive to disk emission.

Our results yield longer disk decay times than those reported by Fedele et al. (2010), who estimate dust dispersion timescales around 3 Myr and mass accretion lifetimes of 2.3 Myr. They argue that accretion ceases earlier than dust dispersal, potentially due to planet formation or migration in the inner disk. Our longer decay times suggest that dust (and, by extension, the disk material) could persist well beyond the cessation of accretion, allowing for an extended window of planet formation. For comparison, Delfini et al. (2025) derive accretion decay timescales in Sco-Cen that are consistent with the values reported by Fedele et al. (2010).

Another important caveat relates to stellar mass. Numerous studies have shown that disk lifetimes vary with mass, with higher-mass stars typically losing their disks more rapidly than their lower-mass counterparts (Carpenter et al. 2006; Roccatagliata et al. 2011; Luhman & Mamajek 2012; Yasui et al. 2014; Ribas et al. 2015). Our analysis does not explicitly stratify by mass but instead marginalizes over the full stellar mass range sampled by the SigMA catalog and Gaia selection function at roughly 100–200 pc distances – spanning from the hydrogen-burning limit up to ∼8 M_⊙. However, due to the steep nature of the initial mass function (IMF), our sample is dominated by low-mass stars, even though we are sensitive to higher masses as well. The resulting disk decay time thus reflects a weighted average over this low- to intermediate-mass regime and could differ in stellar populations with a different mass distribution or in environments with distinct star formation conditions. This might explain the discrepancy to previous studies, as noted by Pfalzner et al. (2022), that might have underestimated disk lifetimes due to sample biases. Future work incorporating spectral types or mass bins will be valuable for isolating mass-dependent trends in disk evolution.

Finally, Michel et al. (2021) find longer disk lifetimes (up to 8 Myr) in low-UV environments, attributing this to reduced external photoevaporation. Although our estimate for Sco-Cen is slightly lower, it remains broadly consistent with their result. Notably, our sample includes ten stellar subpopulations within Upper Scorpius spanning an age range of 3–19 Myrs, which Michel et al. (2021) treat as a single population.

6 Conclusions

This study introduces the first homogeneous sequence of disk fractions as a function of age for an individual star-forming region. We present the initial sequence of the disk fraction in relation to age for one coherent star-forming region. We compiled disk fraction measurements for 33 stellar clusters in the nearby Sco-Cen OB association and analyzed their functional relationship with age. We find that this relation is well described by an exponential decay function over the age range of approximately 3–21 Myr.

To quantify the characteristic disk lifetime, we developed a probabilistic fitting framework that incorporates observational uncertainties and explores a range of prior assumptions and disk selection methods. Our analysis yields a robust estimate of the disk decay time: 5.8 ± 0.3 Myr. This result remains consistent across different model variants, including those with uninformative priors, and when using alternative disk selection methods. When considering all systematic effects, the disk decay time lies between 4 and 7 Myr. Furthermore, our findings are stable even when incorporating disk fraction measurements from the literature or excluding potential binary systems based on Gaia RUWE values. Taken together, these results support a longer disk lifetime than many previous studies, with important implications for the timescales available for planet formation in low-extinction environments such as Sco-Cen.

Data availability

Tables C.2 and C.3 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/707/A216

Acknowledgements

We thank the anonymous referee for their insightful comments that helped to improve the manuscript. Co-funded by the European Union (ERC, ISM-FLOW, 101055318). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the Euro-pean Union or the European Research Council. S.R. acknowledges funding by the Federal Ministry Republic of Austria for Climate Action, Environment, Energy, Mobility, Innovation and Technology (BMK, https://www.bmk.gv.at/) and the Austrian Research Promotion Agency (FFG, https://www.ffg.at/) under project number FO999892674. S. Ratzenböck performed this work as an SAO postdoctoral fellow, and we acknowledge the Smithsonian Institution for their support. J.G. acknowledges funding from the European Union, the Central Bohemian Region, and the Czech Academy of Sciences, as part of the MERIT fellowship (MSCA-COFUND Horizon Europe, Grant agreement 101081195). This work has made use of (1) data from the European Space Agency (ESA) mission Gaia https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement, (2) the Wide-field Infrared Survey Explorer, All-WISE makes use of data from WISE, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration and (3) data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The work has used Python, https://www.python.org/; arviz Kumar et al. (2019), Astropy:http://www.astropy.org, a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022), corner (Foreman-Mackey 2016), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), pandas (pandas development team 2024), pymc (Oriol et al. 2023), scikit-learn (Pedregosa et al. 2011), SciPy (Virtanen et al. 2020), seaborn (Waskom 2021) and Uncertainties: a Python package for calculations with uncertainties, Eric O. LEBIGOT, http://pythonhosted.org/uncertainties/. This research has made use of TopCat (Taylor 2005). This research has made use of the SVO Filter Profile Service “Carlos Rodrigo”, funded by MCIN/AEI/10.13039/501100011033/ through grant PID2023-146210NB-I00.

References

Appendix A Selecting disk candidates

In the following subsections, we outline three different approaches to select young stellar object (YSO) candidates, hence sources with IR excess, to determine the disk fraction of the stellar clusters in Sco-Cen. First and foremost, we use a combination of IR color-color diagrams (CCDs) using WISE and 2MASS photometry (Appendix A.1), second we use the spectral energy distribution (SED) using similar band combinations (Appendix A.2), and, finally, we use a disk selection from the literature (Luhman 2022), to compare to our results using our own disk selection (Appendix A.4).

Appendix A.1 CCD Selection

First, we select disk candidates by using a combination of four CCDs that are composed of 2MASS and WISE photometry, referred to as the CCD selection. This allows us to separate sources with and without IR excess, which is indicative for sources with disks and without disks. We are using a well defined Gaia selected sample of nearby stellar clusters where we can assume that this is free of extragalactic contamination. Therefore, we do not need specific criteria to remove such contaminating sources, which tend to have similar colors as YSOs. We only apply basic quality criteria to identify inferior photometry, as given in the Equs. (A.1) to (A.6). The cross-match of the 12 873 Gaia selected Sco-Cen members with AllWISE and 2MASS IR data yields 10 087 matched sources (78.4 %). After applying the quality criteria, which are outlined below, we are left with 10 016 sources (77.8 % of the whole sample or 99.3 % of the IR sample). We use these 10 016 sources to perform our YSO selection steps and to derive the disk statistics. As a last step in the CCD selection, we investigated several additional color-magnitude diagrams (CMDs) to check for bright sources with clear IR excess that might have been missed by some of our CCD selection cuts.

In the overall selection, we used all four bands from WISE (W1 to W4) and the H and K_S band from 2MASS. The quality criteria for each band are as follows: $$\text{S}/\text{N}(W1)>10,{\chi }^2(W1)<20,\text{e}(W1)<0.2,$$(A.1) $$\text{S}/\text{N}(W2)>10,{\chi }^2(W2)<20,\text{e}(W2)<0.2,$$(A.2) $$\text{S}/\text{N}(W3)>5,{\chi }^2(W3)<20,\text{e}(W3)<0.2,$$(A.3) $$\text{S}/\text{N}(W4)>5,{\chi }^2(W4)<20,\text{e}(W4)<0.2,$$(A.4) $$\text{e}\left(K_S\right)<0.1,$$(A.5) $$\text{e}(H)<\mathrm{0.1.}$$(A.6)

Equs. (A.1) to (A.4) are the quality cuts for the WISE bands and S/N stands for the signal-to-noise ratio, the second cut is made in the reduced χ²-value and the third is the photometric uncertainty in mag. For the measurements in 2MASS only a cut in the photometric uncertainty is made. The quality cuts are applied to the 10 087 sources on the bands that are used for a given CCD selection, as shown in Fig. A.1. Hence, if a source passes the quality cuts used for the respective bands in a CCD, then it is included in this CCD selection (e.g., the WISE W123 selection only includes the cuts of Equs. (A.1) to (A.3)), while one source can appear in multiple CCD selection. A source is considered disk-bearing if IR excess is detected in one of the four CCDs (or CMDs), independently of the results in the other CCDs. For the four CCD selections applied in this work, we used the following band combinations:

• W123 selection: W1 − W2 versus W2 − W3,

• W124 selection: W1 − W2 versus W2 − W4,

• HKW2 selection: H − K_S versus K_S − W2,

• HKW3 selection: H − K_S versus K_S − W3.

The CCDs and the selection borders are shown in Fig. A.1. Additionally, several CMDs are checked, while one CMD is always used in combination with one of the above CCDs:

• W13 CMD: W3 versus W1 − W3 (for the W123 selection),

• W14 CMD: W4 versus W1 − W4 (for the W124 selection),

• HW2 CMD: W2 versus H − W2 (for the HKW2 selection),

• HW3 CMD: W3 versus H − W3 (for the HKW3 selection).

The different band combinations should ensure that sources, which might have inferior photometry in one band, are still included in our disk selection process. The additional CMD checks should ensure that we are not missing obvious, bright YSO candidates, that might have been excluded by the selection borders in the CCDs. We find that only the W14 CMD, concerning the W124 selection, delivers additional disk candidates (hence, the other three CMDs are not further explained in detail below). More details on each of the four disk selections are given in the following subsections. The resulting number statistics of the four CCD selections are given in Table A.1.

Appendix A.1.1 W123 selection

The first CCD is the W123 CCD as it gives the largest number of sources with a clear separation when using MIR colors. This band combination is also often used in the literature (see e.g. Koenig & Leisawitz 2014; Großschedl et al. 2019, and also Spitzer selections, Evans et al. 2009). We apply quality criteria to the IR photometry using the quality cuts Eq. (A.1) to (A.3). As the sources are preselected using Gaia, we do not need to apply special cleaning steps to avoid extragalactic contamination, and we follow a simpler and similar approach as chosen in Großschedl et al. (2021). Among the sources, which pass the quality cuts, we select sources as disk-bearing if they meet the following requirement in the CCD of W1 − W2 versus W2 − W3: $$W1-W2>-4\cdot (W2-W3)+\mathrm{4.2.}$$(A.7)

The separation is based on information from known YSO samples and their colors. The selection slope of the border is parallel to the extinction vector, informed by the locations of known YSOs and by YSO selection from the literature (see for example Koenig & Leisawitz 2014; Großschedl et al. 2019). The extinction law for the IR bands is taken from Meingast et al. (2016) and Großschedl et al. (2019).

Appendix A.1.2 W124 selection

The W124 CCD is used to select sources with strong IR excess in W4. A similar approach is presented in Koenig & Leisawitz (2014), where they use this band combination to select sources that are potential transition disks (see also Teixeira et al. 2012). The quality criteria are Eq. (A.1), (A.2) and (A.4). To inform our decision for the selection borders in W124 we use the experience from the literature (e.g. Teixeira et al. 2012; Koenig & Leisawitz 2014; Großschedl et al. 2019). Moreover, we use the already selected disk sources from the W123 selection (with additional W4 quality criteria, see Eq. (A.4)), and overplot them in the W124 color space to identify the typical colors of disk sources. Based on this and the literature, two conditions for disk-bearing sources are constructed. $$\begin{array}{l}W1-W2<2\cdot (W2-W4)-3.5,\\ W1-W2>\mathrm{0.1.}\end{array}$$(A.8)

We include this additional horizontal cut in W1 − W2 as we believe there is still some contamination in the W4 band after the quality criteria, which would appear as excess in W4, but not in W1 − W2. Therefore, we cut at the bottom of the CCD to avoid this contamination.

The CMD of W4 versus W1 − W4 shows bright sources with clear IR excess, which are classified as disk-less sources in the W124 CCD. Therefore, we add sources as disk-bearing if they fulfill: $$\begin{array}{l}W4<4,\\ W1-W4>4.\end{array}$$(A.9)

With the W14 CMD, we added four additional disk candidates.

Appendix A.1.3 HKW2 selection

The HKW2 CCD is used additionally to select sources with inferior W3 or W4 photometry, or where these two longer wavelength bands are influenced by crowding. We apply the quality cuts from Equs. (A.2), (A.5) and (A.6), and sources fulfilling the following selection cut are selected as disk-bearing candidates: $$H-K_S<\frac{0.55}{0.45}\cdot \left(K_S-W2\right)-\mathrm{0.75.}$$(A.10)

Again the slope is parallel to the extinction vector. The separation is found by plotting the selection of disk-bearing and disk-less sources from the W123 CCD in the HKW2 CCD, while only sources that pass the additional HKW2 quality criteria are used from the W123 selection for this comparison. The separation is again informed by locations of YSOs from the literature (e.g. Teixeira et al. 2012).

Appendix A.1.4 HKW3 selection

Finally, for the HKW3 selection, we apply the quality criteria from Equs. (A.3), (A.5) and (A.6). Sources, which fulfill the following selection criteria, are added as disk-bearing candidates: $$H-K_S<\frac{0.55}{0.39}\cdot \left(K_S-W3\right)-\mathrm{2.3.}$$(A.11)

Again, the selection border is chosen parallel to the extinction vector and by using the information of the previously selected disk candidates or experience from the literature (see also the preceding sections).

Appendix A.2 SED Selection

We present an alternative selection approach by using the not extinction corrected MIR spectral energy distribution (SED) instead of the CCDs, to have an independent list of YSO candidates to be used to estimate the disk fraction per stellar cluster. With this we will compare the selections among each other and use both selections in further analysis to see what influence the used selection method might have on our final results. We call this selection the SED selection.

The magnitude values are taken with their corresponding uncertainties and converted into flux density in terms of wavelength F_λ. For the flux zero-points and central wavelengths of the bands, we use the SVO Filter Profile Service (Rodrigo et al. 2012; Rodrigo & Solano 2020; Rodrigo et al. 2024). The selection is done using the spectral index (Lada 1987). $$\alpha =\frac{\text{d}\log \left(\lambda F_{\lambda }\right)}{\text{d}\log \lambda }.$$(A.12)

To determine α, a linear fit of the logarithm of the flux density times the wavelength versus the logarithm of the wavelength needs to be performed. The slope of this function is taken as α and then used for selecting disk candidates.

In the literature (Lada et al. 2006; Evans et al. 2009; Teixeira et al. 2012; Dunham et al. 2015; Großschedl et al. 2019), the Spitzer IRAC bands or the 2MASS K_S band together with the IRAC bands are frequently used to define the YSO classes using the spectral index. We strive to use similar band combinations, and we use the range from K_S to W3 as our default band range to calculate the spectral index α, ideally including all bands in between (W1 and W2). We require that the end points of the band range pass the quality criteria. The bands in the middle do not have to pass the quality criteria, but a measurement has to exist. If a source does not pass the quality criteria of the K_S or W3 (hence, if K_S or W3 is missing), then we use another band combination, given in Table A.2. This could lead to some inconsistencies in the YSO selection, however, it is the best approach to gain additional candidates, mostly for the cases where W3 is missing.

Fig. A.1 Selection of disk candidates obtained from a combination of four different IR CCDs. The upper left panel shows the W123 CCD selection (W1 − W2 versus W2 − W3), the upper right panel the W124 CCD selection (W1 − W2 versus W2 − W4), the lower left panel the HKW2 CCD selection (H − KS versus KS − W2), and the lower right the HKW3 CCD selection (H − KS versus KS − W3). In each plot, the blue circles show sources selected as disk-bearing and the orange crosses are disk-less. The borders are the dashed, black lines in each diagram. In each plot, the extinction vector is shown, scaled to a value of AKS given in the plot above the extinction vector.

Finally, the SED disk selection is defined as follows: $$\alpha \ge -2.$$(A.13)

The remaining sources that are passing the quality criteria of one of the used band combinations in Table A.2 are defined as diskless. The cut at −2 is based on the location of known disk sources in the SED slope distribution in Fig. A.2. Moreover, we see in the distribution that there is a minimum around −1.8. We do not use this value for our cut, but we decide on the slightly less conservative −2 threshold. This is also supported by the literature. Lada et al. (2006) find that pre-main-sequence stars without a circumstellar disk (hence, the unobscured photospheres) should have SED slopes of about −2.6 (when using the IRAC bands). We pick −2 to account for systematic uncertainties and scatter due to measurement uncertainties.

Figure A.2 shows the histogram of the SED slopes (α) of all sources for which an SED fit is possible in one of the band combinations, given in Table A.2. One can see two peaks, roughly representing the disk-bearing and disk-less sources, while the disk-less sources dominate in the whole Sco-Cen OB association, which includes clusters up to about 21 Myr. The statistics of the SED fitting for the different combinations of bands and the total numbers of the SED selection are presented in Table A.2.

Fig. A.2 Histogram of the spectral index α for all sources obtained through SED fitting. The orange bars indicate disk-less sources, the blue one disk-bearing. The selection border at −2 is shown as black, dashed line.

Fig. A.3 Confusion Matrix for the selection of disk-bearing and diskless sources using the CCD selection (y-axis) and the SED selection (x-axis). We give the percentages in each cell relative to the number common in both selection methods (10 008 sources). The cells are additionally color-coded by the percentage value.

Appendix A.3 Comparing the CCD and SED selections

To assess the two presented selection methods (CCD and SED), we construct a confusion matrix between the two, which is shown in Fig. A.3. Among the 10 016 sources in the CCD selection and 10 021 sources in the SED selection, 10 008 sources appear in both. We divide the total number in each cell of the confusion matrix by the absolute number of sources (10 008 sources), to give the percentages of the confusion matrix. The main diagonal represents the same selection in both selection methods while the off-diagonal entries show different selections. In Table A.1 and A.2 we already see that the percentage of disk-bearing and disk-less sources is almost the same. We also see this in the confusion matrix, but there is a slight degree of confusion between disk-bearing and disk-less sources. This confusion is very similar, which is why the fraction of disk-bearing and disk-less sources is the same for the CCD and SED selection method.

Fig. A.4 The absolute difference in disk fraction f for the selection using the CCD versus the SED or the Luhman selections, plotted against stellar cluster age. The light red data points are the 33 stellar cluster of the CCD/SED selection and the light green squares the 25 stellar clusters that are included in the Luhman selection. The histogram on the right along the y-axis shows the distributions of the difference for the two comparisons.

Computing the disk fraction for each stellar cluster using the CCD and SED selections yields slightly different values per stellar cluster. This is because the similar number of disk-bearing and disk-less sources does not have to be distributed the same way among the stellar clusters. To study the effect, we plot the absolute difference of disk fraction using the CCD selection minus the SED selection versus the age, which is shown in Fig. A.4 (light-red dots and histogram). We see that the data points scatter randomly around the zero-line (grey dashed line). We perform a linear fit and calculate the R²-value to assess if the relation is better described by a linear or a constant relation. The R²-value of the linear fit yields 0.008, so the absolute difference in disk fraction and age is best described through a constant relation. Calculating the mean and the standard deviation of the difference in disk fraction between CCD and SED selection yields 0.005 ± 0.019, so the difference is negligible.

Appendix A.4 Comparison to Luhman (2022)

Additionally, we compare our CCD selection to results from Luhman (2022). Luhman (2022) classifies the sources in Sco-Cen based on IR excess into different disk types. We collapse this classification into a broad disk-bearing and disk-less separation and refer to this as the Luhman selection. We cross-match their table to our data of 12 873 sources based on the Gaia and WISE identifiers and recover 8137 sources (63.2 %). We remove empty classifications and sources classified as Be stars. The following of their classes are then considered as disk-bearing: full, evolved, transitional, ev or trans. Also sources with uncertain classification (indicated by a ? or edge-on?) are included. Sources classified as debris/ev trans, debris or III (with and without ?) are considered as disk-less. We obtain 8137 sources, where 920 (11.3 %) are disk-bearing and 7217 (88.7 %) are disk-less. We see a bit lower percentage of disk-bearing sources in comparison to the CCD and SED selection (see Table A.1 and A.2).

We compare the Luhman selection similarly to the CCD selection and we construct a confusion matrix, as shown in Fig. A.5. 7889 sources appear in both the CCD and Luhman selection. Again we see a slight confusion in the off-diagonal cells, which is comparable to the confusion matrix of CCD and SED selection in Fig. A.3. The percentage of disk-bearing sources is overall smaller for this sample, which could be due to the missing stellar clusters, which are usually younger clusters.

Fig. A.5 Confusion Matrix for the selection in disk-bearing and diskless sources using CCDs (y-axis) and in comparison to Luhman (2022) (x-axis). We give the percentages in each cell relative to the number common in both selection methods (7889 sources). The cells are additionally color-coded by the percentage value.

We again compute the disk fraction of each of the 25 available stellar clusters based on the Luhman selection and compare this to the results from the CCD selection, which is applied on the full SigMA cluster sample. The absolute difference versus stellar cluster age is shown in Fig. A.4 (light green squares and histogram). Again, we perform a linear fit to see if the relation is better described by a linear relation or a constant relation. The R²-value yields 0.017, so again it is better described by a constant relation. The mean and standard deviation of the difference in disk fraction is 0.003 ± 0.024, again a negligible variation.

Appendix B Statistical model

Appendix B.1 Mathematical motivation

With the exponential model function defined in Eq. (2), an estimated disk fraction $\widehat{f}\left(t_i\right)$ can be calculated for a stellar cluster i at time t_i, where i ∈ {1, …, 33}. Assuming the number of sources N_i in stellar cluster i is fixed, the estimated number of disks ${\widehat{k}}_i$ in stellar cluster i becomes $${\widehat{k}}_i=N_i\cdot \widehat{f}\left(t_i\right)=N_i\cdot f_0\cdot e^{-\frac{t_i-t_0}{\tau }}.$$(B.1)

The observed disk fraction f_i of stellar cluster i is affected by random additive noise e_i. Thus, the true disk fraction y_i is related with the observed value by $$f_i=y_i+e_i.$$(B.2)

In the following we aim to motivate the model for the noise e_i. The noise model has to reflect the underlying counting process, namely Bernoulli, where the disk fraction of a given stellar cluster governs the successfully identified disks k_i out of a sample of N_i sources. Hence, for a given disk fraction f_i the number of detected disks k_i out of a total of fixed N_i sources is binomially distributed. $$p\left(k_i\mid f_i,N_i\right)=Binomial\left(N_i,f_i\right).$$(B.3)

Instead of studying p(k_i | f_i, N_i), we have data of k_i and N_i and aim to investigate p( f_i | k_i, N_i) which quantifies the disk fraction uncertainty e_i. We invert this statement using Bayes formula, $$p\left(f_i\mid k_i,N_i\right)=\frac{1}{Z}\cdot p\left(k_i\mid f_i,N_i\right)\cdot p\left(f_i\mid N_i\right).$$(B.4) $$Z={\displaystyle \int d}{f}_{i}p\left(k_i\mid f_i,N_i\right)p\left(f_i\mid N_i\right).$$(B.5)

Here, Z is a normalization factor that guarantees p( f_i | k_i, N_i) is a proper probability distribution that we can ignore. We model p( f_i | N_i) to be independent of N_i and as being a Beta distribution with parameters α = 1 and β = 1, producing a uniform distribution over the interval (0, 1). The Beta distribution is a conjugate distribution of the binomial distribution leading to an analytically tractable compound distribution in which the nor-malization integral needs not to be manually evaluated. The compound distribution p( f_i | k_i, N_i) turns out to be the following Beta distribution. $$p\left(f_i\mid k_i,N_i\right)=p\left(e_i\right)=Beta\left(\alpha =k_i+1,\beta =N_i-k_i+1\right).$$(B.6)

We can write down the likelihood of seeing a given disk fraction conditioned on the physical parameters, $$\begin{array}{l}p\left(f_i\mid \tau ,t_0,f_0\right)=\hfill \\ Beta\left(\alpha ={\widehat{k}}_i+1,\beta =N_i-{\widehat{k}}_i+1\right)=\hfill \\ Beta\left(\alpha =N_i\cdot f_0\cdot e^{-\frac{t_i-t_0}{\tau }}+1,\beta =N_i\cdot \left(1-f_0\cdot e^{-\frac{t_i-t_0}{\tau }}\right)+1\right),\hfill \end{array}$$(B.7) where the estimated number of disks from Eq. (B.1) is used. Assuming independence between data points, the probability distribution function (PDF) for the full data set L_S is the product of single observations, $$\begin{array}{l}{ℒ}_S=p\left({\left\{f_i\right\}}_{i=1}^{33}\mid \tau ,t_0,f_0\right)=\hfill \\ ={\displaystyle \prod _{i=1}^{33}p}\left(f_i\mid \tau ,t_0,f_0\right).\hfill \end{array}$$(B.8)

In addition to modeling the signal behavior, we also aim to mitigate the effect of outliers and deal with unaccounted noise sources in the observed disk fraction. The contribution of outlier sources is aimed to be modeled directly as a separate component in the likelihood using a mixture model approach. From a visual inspection of the disk fraction over age (see Fig. 2), we find a couple of outliers which appear to lie below the common exponential trend line, meaning a lower disk fraction. As the number of disks in these clusters is very low, their upper uncertainty is usually also low (in contrast to points above the common exponential trend line). Assuming outliers are age independent and group towards small disk fractions, we model their contribution to the likelihood ℒ_B via a Half-Normal distribution with a set mean of µ = 0 and variable standard deviation, determined during parameter inference, $${ℒ}_B=\mathcal{N}(0,\sigma ).$$(B.9)

Thus, this component aims to avoid biasing the exponential fit towards smaller disk fractions. These two probability distributions add up to the combined likelihood ℒ as $$ℒ=w_1{ℒ}_S+w_2{ℒ}_B,$$(B.10) where w₁ and w₂ are the weights that fulfill w₁ + w₂ = 1.

Appendix B.2 Implementation

Using Bayes’ theorem the posterior PDF of the model parameters in Eq. (2) can be calculated using the combined likelihood, $$p\left(\tau ,t_0,f_0\mid {\left\{f_i\right\}}_{i=1}^{33}\right)\propto ℒ\cdot p\left(\tau ,t_0,f_0\right).$$(B.11) p(τ, t₀, f₀) is the prior assumption of the parameters in form of a PDF. Prior probability distributions are chosen for the parameters of the exponential fit.

The prior for the standard deviation σ of the background distribution in Eq. (B.9) is a truncated Normal distribution between 0.01 and 0.5, with 0.1 being the value for mean and its standard deviation. We tested the fitting approach for the fiducial model without incorporating the background component into a joint mixture model (i.e., using only the likelihood given by ℒ_S). This changes the posterior median of the decay time to ${6.16}_{-0.20}^{+0.18}$ Myr.

We build and fit the model using the python library pymc (Oriol et al. 2023). In the case of the stellar cluster having zero disk-bearing sources (and with that a disk fraction of zero), we change its value to 10⁻⁷ to guarantee numerical stability. The sampling is done using NUTS (No U-Turn Sampler) (Hoffman et al. 2014). During sampling divergences can occur, this is mitigated by tuning the target accept probability so that the model returns 0 divergences, which is found at a value of 0.99. To assess the quality of the resulting samples, we compute the rank normalized $\widehat{R}$ value for all fitted parameters and check if it is below 1.01 (Vehtari et al. 2021). This is fulfilled for all model results presented here.

For each parameter, we compute its posterior median and determine its 1-σ confidence interval by computing the 68 % highest density interval (HDI).

The described process of modeling the counting of disks using a Bernoulli distribution can be used to obtain uncertainties for the disk fractions. With this, we obtain 1-σ confidence intervals for the disk fractions, shown in Fig. 2.

Appendix B.3 Prior choices for the parameters

We fit a total of three different models to the data. These differ in the modeling (prior PDFs or functional relation). Two of the fitting results differ in the selection approach used and two results are subsets of the CCD selection data set.

Our reference statistical model, the fiducial model, aims to include physical knowledge of the parameters in its priors. The prior PDF for the decay τ is a Half-Cauchy distribution with scale parameter β = 5. The value is chosen to initiate a flat prior and the initial guess is set to 5, which is determined through least-squares fitting of Eq. (2). For the intercept f₀ a Beta distribution is used as the value is confined between 0 and 1. The shape parameters are chosen as α = 5 and β = 1 to put more weight on values closer to 1. The prior for the shift t₀ is a truncated Normal distribution with mean 0.1 and standard deviation 10. As the shift should move the function towards the data, the lower border is 0 and the upper 10. The borders are constructed by shifting the value found through the least-squares fit. The mean and standard deviation are chosen in a way that the prior PDF is flat.

To study the influence of the physical assumptions made on the priors, we also fit a data-driven model where no constraint is put on the model priors and the optimal values are determined by the data. All parameters here are modeled using truncated Normal distributions. The prior of the decay τ has a mean of 5, due to result from least-squares fitting, and standard deviation of 10 to initiate a wide distribution. The borders are 0 and ∞ as the data indicates exponential decay and not growth. The intercept f₀ prior has as mean 1 and standard deviation 0.5. This is chosen to study the behavior around the physical border of the intercept value. The lower border is 0 and the upper 10 to allow the exploration of larger intercept values. The shift t₀ is centered on 0 with a standard deviation of 5. Here we want to study if the data would prefer a shift of the function to the left (represented by negative values of t₀). The borders are −5 and 5.

The shift influences the decay and we want to see its influence on the parameters. This is why we fit the data-driven model again, but the function in the model is changed from Eq. (2) to $$f(t)=f_0\cdot e^{-\frac{t}{\tau }},$$(B.12) by setting t₀ to 0. This model is called the no-shift model.

The selection is also done by SED fitting (SED selection, see Appendix A). Their results are converted to disk fractions for the stellar clusters and the fiducial model is fitted to the data.

In Appendix A, we compare our selection to the one of Luhman (2022) (Luhman selection). For their selection, we compute the disk fraction for each stellar cluster and fit the fiducial model to them.

As a final test we try to deduce the effect of binaries or multiples on the disk fraction (see for example Kraus et al. 2012). To this end we use the Gaia RUWE parameter, while higher RUWE values (RUWE ≥ 1.4) indicate that a source might be a stellar binary or multiple candidate (Penoyre et al. 2022a,b; Castro-Ginard et al. 2024). Hence, we apply a quality selection on RUWE using the CCD selection. The first subsample is RUWE < 1.4 and the second is RUWE ≥ 1.4. We follow the choice of Ratzenböck et al. (2023b) for the threshold value of 1.4 to define this separation consistently, where they use it to exclude binary candidates when inferring cluster ages through isochrones. We again determine the disk fractions individually for the two subsets, and then we apply the fiducial model on the disk fraction versus cluster age distribution. The RUWE ≥ 1.4 subset is mainly used for comparison purposes, to better understand the potential influence of binaries on the evaluation of disk decay times.

Appendix C Auxiliary Figures and Tables

In Fig. C.1 to C.3 we show corner plots for the fiducial, the data-driven and the no-shift model. In Fig. C.4 we show the same as Fig. 2 using different samples of stellar cluster ages and fidu-cial+CCD results for the different plots. In Table C.1 we show the stellar cluster statistics, number of disk-bearing and disk-less sources and disk fraction for the different selection methods.

Table C.2 provides the column description for the table of stellar members, based on the stellar member catalog from Ratzenböck et al. (2023a). Table C.3 details the description for stellar clusters in Sco-Cen as used in this work.

Fig. C.1 Corner plot of the sampled parameters as determined with the fiducial model (see Appendix B). The left upper panel shows the distribution of the decay parameter τ in Myr. In the second row, the first panel shows the correlation of τ with the shift parameter t0, and the second panel shows the distribution of t0. The last row shows first the correlation of τ versus the intercept parameter f0, then the correlation of t0 versus f0, and finally the distribution of the intercept parameter. The red line shows the median of each sampled parameter.

Fig. C.2 Corner plot of the sampled parameters as determined with the data-driven model (see Appendix B). The structure and color are the same as in Fig. C.1.

Fig. C.3 Corner plot of the sampled parameters as determined with the no-shift model (see Appendix B). The structure and color are the same as in Fig. C.1 while the shift parameter t0 is excluded here, since it is a constant in the no-shift model.

Fig. C.4 Disk fraction versus age for different samples of ages. The black data points again represent the stellar clusters, where the disk fraction is constant among all plots (CCD selection). The top left plot shows ages from PARSEC-G-RP, the top right Baraffe-G-RP, the bottom left Baraffe-BP-RP and the bottom right the ages from Kerr et al. (2021). The orange, dashed line denotes the median of the fiducial+CCD result applied on these ages and the gray area is the 1σ confidence interval. The red solid line gives the fiducial+CCD result from Fig. 2 using the PARSEC-BP-RP ages. The uncertainties in ages are from Ratzenböck et al. (2023b) or Kerr et al. (2021) and the disk fraction uncertainties are the same as in Fig. 2.

All Tables

All Figures

	Fig. 1 Plot of galactic longitude versus galactic latitude showing the region of Sco-Cen. The sources are shown based on their selection, blue circles for disk-bearing, and orange crosses for disk-less sources.
In the text

Fig. 2 Disk fraction versus age. The black data points represent the stellar clusters, on the left for the CCD selection, in the middle for the SED selection and on the right for the Luhman selection. The Luhman selection only includes 25 out of 33 stellar clusters because the study by Luhman (2022) encompasses a smaller region compared to Ratzenböck et al. (2023a). Shown in green-blue in the right plot are the missing stellar clusters (from the CCD selection). In the left plot, the red line represents the median of the fiducial+CCD result, which is also shown in the middle and right plot with reduced transparency. The orange dash-dotted line in the middle plot gives the median of the fiducial+SED result. The blue dotted line is the median of the fiducial+Luhman result. The gray area shows the 1σ confidence interval for each median result. The uncertainties in ages are from Ratzenböck et al. (2023b) and the disk fraction uncertainties are derived from the underlying counting process which can be modeled as a simple Bernoulli trial.

In the text

Fig. 3 Violin plot showing the posterior probability density functions (PDFs) of the disk decay time parameter τ. The x-axis denotes the decay time in Myr, while the y-axis labels each model and disk selection method. The number N displayed next to each distribution indicates the number of sources used in the corresponding fit; this number remains constant for the top three rows, where only the model configuration varies. Posterior distributions are grouped and color-coded: red indicates variations in prior assumptions (with the fiducial model highlighted at the top), green represents the supplementary disk selection methods (SED, Luhman) and blue corresponds to subsamples based on stellar multiplicity (single vs. candidate binaries). See Table 1 for the corresponding numerical results.

In the text

Fig. A.1 Selection of disk candidates obtained from a combination of four different IR CCDs. The upper left panel shows the W123 CCD selection (W1 − W2 versus W2 − W3), the upper right panel the W124 CCD selection (W1 − W2 versus W2 − W4), the lower left panel the HKW2 CCD selection (H − KS versus KS − W2), and the lower right the HKW3 CCD selection (H − KS versus KS − W3). In each plot, the blue circles show sources selected as disk-bearing and the orange crosses are disk-less. The borders are the dashed, black lines in each diagram. In each plot, the extinction vector is shown, scaled to a value of AKS given in the plot above the extinction vector.

In the text

	Fig. A.2 Histogram of the spectral index α for all sources obtained through SED fitting. The orange bars indicate disk-less sources, the blue one disk-bearing. The selection border at −2 is shown as black, dashed line.
In the text

	Fig. A.3 Confusion Matrix for the selection of disk-bearing and diskless sources using the CCD selection (y-axis) and the SED selection (x-axis). We give the percentages in each cell relative to the number common in both selection methods (10 008 sources). The cells are additionally color-coded by the percentage value.
In the text

Fig. A.4 The absolute difference in disk fraction f for the selection using the CCD versus the SED or the Luhman selections, plotted against stellar cluster age. The light red data points are the 33 stellar cluster of the CCD/SED selection and the light green squares the 25 stellar clusters that are included in the Luhman selection. The histogram on the right along the y-axis shows the distributions of the difference for the two comparisons.

In the text

	Fig. A.5 Confusion Matrix for the selection in disk-bearing and diskless sources using CCDs (y-axis) and in comparison to Luhman (2022) (x-axis). We give the percentages in each cell relative to the number common in both selection methods (7889 sources). The cells are additionally color-coded by the percentage value.
In the text

Fig. C.1 Corner plot of the sampled parameters as determined with the fiducial model (see Appendix B). The left upper panel shows the distribution of the decay parameter τ in Myr. In the second row, the first panel shows the correlation of τ with the shift parameter t0, and the second panel shows the distribution of t0. The last row shows first the correlation of τ versus the intercept parameter f0, then the correlation of t0 versus f0, and finally the distribution of the intercept parameter. The red line shows the median of each sampled parameter.

In the text

	Fig. C.2 Corner plot of the sampled parameters as determined with the data-driven model (see Appendix B). The structure and color are the same as in Fig. C.1.
In the text

	Fig. C.3 Corner plot of the sampled parameters as determined with the no-shift model (see Appendix B). The structure and color are the same as in Fig. C.1 while the shift parameter t0 is excluded here, since it is a constant in the no-shift model.
In the text

Fig. C.4 Disk fraction versus age for different samples of ages. The black data points again represent the stellar clusters, where the disk fraction is constant among all plots (CCD selection). The top left plot shows ages from PARSEC-G-RP, the top right Baraffe-G-RP, the bottom left Baraffe-BP-RP and the bottom right the ages from Kerr et al. (2021). The orange, dashed line denotes the median of the fiducial+CCD result applied on these ages and the gray area is the 1σ confidence interval. The red solid line gives the fiducial+CCD result from Fig. 2 using the PARSEC-BP-RP ages. The uncertainties in ages are from Ratzenböck et al. (2023b) or Kerr et al. (2021) and the disk fraction uncertainties are the same as in Fig. 2.

In the text