HST Survey of the Orion Nebula Cluster in ACS/Visible and WFC3/IR Bands. IV. A Bayesian Multiwavelength Study of Stellar Parameters in the Orion Nebula Cluster

We have performed a comprehensive study of the Orion Nebula Cluster (ONC) combining the photometric data obtained by the two Hubble Space Telescope Treasury programs that targeted this region. To consistently analyze the rich data set obtained in a wide variety of filters, we adopted a Bayesian approach to fit the spectral energy distribution of the sources, deriving mass, age, extinction, distance, and accretion for each source in the region. The three-dimensional study of mass distribution for bona fide cluster members shows that mass segregation in the ONC extends to subsolar masses, while the age distribution strongly supports the idea that star formation in the ONC is best described by a major episode of star formation that happened ∼1 Myr ago. For masses ≳0.1 M ⊙, our derived empirical initial mass function (IMF) is in good agreement with a Chabrier system IMF. Both the accretion luminosity (L acc) and mass accretion rates ( Ṁacc ) are best described by broken power-law relations. This suggests that for the majority of young circumstellar disks in this cluster the excess emission may be dominated by X-ray-driven photoevaporation by the central star rather than external photoevaporation. If this is the case, the slopes of the power-law relations may be largely determined by the initial conditions set at the onset of the star formation process, which may be quite similar between regions that eventually form clusters of different sizes.


INTRODUCTION
The Orion Nebula Cluster (ONC) represents one of the richest (∼ 2000 members in the inner 2 pc radius; McBride & Kounkel 2019) and youngest (∼ 2 Myr; Jeffries et al. 2011;Reggiani et al. 2011;Jerabkova et al. 2019) star-forming regions (SFRs) within 2 kpc (∼ 402 pc; Kuhn et al. 2019) from the sun (Lada & Lada 2003;Portegies Zwart et al. 2010).Due to its proximity and modest foreground extinction (A V ∼ 1 Scandariato et al. 2011), the ONC has been thoroughly studied at visible and near-infrared wavelengths and is therefore re-garded as an ideal laboratory where to investigate, down to very low mass objects, critical aspects of star and planetary formation, such as the initial mass function, mass segregation and early dynamical evolution, radiative feedback and protoplanetary disk photoevaporation and evolution.
In this paper, we combine for the first time the data coming from the two HST Treasury programs (GO-10246 and GO-13826, P.I.M. Robberto) that targeted the ONC.Taking advantage of the extremely accurate photometry provided by the HST in a wide selection of filters, we adopt a Bayesian approach with Markov Chain Monte Carlo (MCMC) strategy to estimate the main properties of each star in the cluster (e.g.mass, extinction and age, distance and accretion) according to the BT-Settl family of isochrones.A Bayesian approach Table 1.ACS visits strategy and photometric system from Robberto et al. (2013).The zero points for each filter are obtained from the STScI ACS Zeropoints Calculator.
allows incorporating as priors other available information e.g. on the source effective temperature, distance, reddening, membership, etc. obtaining a most accurate and consistent estimate of the stellar parameters.
Our first aim is to assess the viability of this technique, considering that our dataset is largely based on archival data, taken at different epochs, and the number of filters available for each source is not uniform across the sample.Concentrating on the most reliable sources, we then perform a deeper analysis of the star formation and evolution history of the Orion Nebula Cluster.
The paper is organized as follows: the observations are presented in §2, including a new ACS photometric catalog based on the most recent instrument calibration data.In §3 we introduce our Bayesian code to perform Spectral Energy Distribution (SED) fitting.In Section §4 we present our final catalog of bonafide cluster members, the derived stellar parameters, and the mass accretion luminosities and rates.In Section §5 we discuss the implication of our findings.Our conclusions are summarized in §6.

DATASETS AND NEW ACS PHOTOMETRY
The HST program GO-10246, executed between October 2004 andMarch 2005, observed the ONC with the F435W, F555W, F658N, F775W and F850LP filters of ACS/WFC, the F336W, F439W, F656N, and F814W of WFPC2, and F110W and F160W of NICMOS (Robberto et al. 2013).The richest dataset is the one provided by ACS, which includes 3399 unique sources, most of them observed two or more times in different HST visits, separated by time intervals ranging from a few hours to several months.Robberto et al. (2013) lists all individual detection, a grand total of 8185 entries each with multi-color observations.In 2015 the HST program GO-13826 returned to the ONC with the F130N and F139M of WFC3 (Robberto et al. 2020).As shown by Figure 1, the areal coverage is similar for the ACS, WFPC2, and WFC3 data, while the NICMOS observations sampled about 1/4 of the field due to the small size of the detector.
In Table 2 we show the number of sources detected in each filter.We will use this dataset as our reference catalog.
Given both the advances in the ACS calibration and the possibility of obtaining more accurate data on close pairs, we reanalyzed the ACS dataset to derive updated aperture photometry.In Table 1 we list the main parameters of the photometric system for the five ACS bandpasses including the most recent assessment of the zero points at the epoch of observations, as provided by the ACS zero points Calculator 1 .Compared to the original zero points implemented by the Robberto et al. (2013), the new values show a difference between −0.002 and 0.019 magnitudes, with the smallest difference in the F555W filter and the highest one for both the F775W and F850LP filters.
Aperture photometry was obtained using the aperture photometry package from the StraKLIP pipeline (Strampelli et al. 2022).This tool allows the detection and subtraction of the signal from sources very close to the primary target to perform accurate photometry on both components.In this paper we shall focus on the main population of isolated stars, leaving the analysis of close pairs to a future paper.
The individual measures obtained in different visits were finally averaged vetting and rejecting those with high uncertainty due to the residual detector defects or the presence of cosmic rays in the immediate vicinity.
In Figure 2 we plot our final averaged magnitudes with their estimated uncertainties for all sources in the catalog in the five ACS filters.Saturation starts at about m 435 = 16, m 555 = 15.75, m 658 = 12.25, m 775 = 15.25 and m 850 = 14.25, corresponding respectively to 0.7, 0.5, 1.2, 0.14 and 0.13 M ⊙ , assuming the BT-Settl 1M yr isochrone at 400 pc, without extinction or accretion.At the other extreme, the 5σ sensitivity limits (σ mag ≃ 0.2) are at m 435 = 23.77,m 555 = 24.07,m 658 = 20.15,m 775 = 23.84 and m 850 = 23.01,roughly corresponding to 0.03, 0.02, 0.04, 0.005 and 0.003M ⊙ , under the same assumptions.limits as a function of extinction and detectable mass.All magnitudes are in the Vega system.Apart from ACS, for the other instruments and passbands we directly adopted the catalogs presented by Robberto et al. (2013Robberto et al. ( , 2020)).
Since most observations were not carried out simultaneously, variability represents a source of uncertainty.According to Herbst et al. (2002), about half of the stars brighter than I ≃ 16 show peak-to-peak variations of ∼ 0.2 mag, or more.While this is typically not the range of magnitudes probed by our deeper HST observations, it suggests that variability randomly affects the quality of the fit.In principle, systematic uncertainties could be reduced using colors combining measures obtained with the same instrument within the same HST orbit, i.e. within about 45 minutes.Testing this approach, we found that in most cases having fewer data points with uncertainties added in quadrature precludes achieving any significant gain compared to the case where the SED fit is performed using single magnitudes.We have  ... 3 is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.Note that in the published machine-readable version of this table, the errors will occupy a different column.

Note-Table
therefore carried out our analysis ignoring variability, at least initially, returning later to the sources that may not pass our tests on the quality of the fit because of spurious photometric data for a specific epoch of observation.
Table 3 shows a sample of our final photometric catalog for all 3399 sources, fully available in electronic format.
The table is organized as follows: the first column reports the ID as provided from the StraKLIP pipeline for cross-identification, while the following columns report the average number of saturated pixels, number of identified bad pixels, magnitude, and uncertainty for each filter.

Bayesian analysis
For our analysis we assume that the main unknown variables that characterize each source are mass, extinction, age, distance, and the accretion parameter (SP acc ), representing the fraction of stellar bolometric luminosity that can be attributed to accretion. .To determine them, we perform Spectral Energy Distribution (SED) fitting comparing the observed photometry to synthetic photometry derived from theoretical models of the BT-Settl family (see Section 3.1.1).In order to obtain a probability distribution for each parameter, our fitting procedure adopts a Bayesian approach with an MCMC algorithm, represented schematically in Figure 3 and based on the following main steps: ): ):

PDJ
)1 ): PDJ )/3 • at the start of the MCMC run, five parent distributions are generated uniformly spanning the parameter space presented in Table 4.The green colored blocks represent fixed ingredients for the simulation while in blue are shown the fitting parameters and in white the key steps of the Bayesian approach.The yellow block shows the stored final output of each iteration.Note that all stars are treated completely independently until the endpoint.
• for each star, a user's defined number of walkers is generated, limited in this work to a maximum of 100.Then the following loop is iterated for each walker (referred to as the "walker loop" in Figure 3): from each distribution discussed above (see also Table 4), a random value is extracted to provide the initial state for the walker in that variable.
-The observed magnitudes of the star under consideration are pulled from the catalog and parsed to the fitting routine as fixed ingredients for the simulation.If one or more filters show sign of saturation, they are discarded from the fit.
-The five parameters randomly extracted from their parent distributions are combined with the models (see Section 3.1.1)to obtain synthetic magnitudes.
-The synthetic and observed magnitudes are iteratively compared to find the set of parameters that better reproduces the observations.This optimization requires: 1) priors on the fitting parameters (see Section 3.1.3),2) a likelihood function to compare observations and simulations (see Section 3.1.2),3) an algorithm for sampling the posterior distribution of the fit parameters (see Section 3.1.4); -The walkers are let to evolve until convergence is reached or a maximum number of steps is executed (2000 in our case).The resulting posterior distribution of the parameters is saved for analysis and the routine moves to the next source; In the following subsections, we detail the main steps of the procedure.

Models
For our analysis, we use the evolutionary models and spectra of the BT-Settl family (AGSS2009 models using the Asplund et al. (2009) solar abundances) 2 , spanning the age range between 1 Myr to 10 Gyr and masses from 0.005 M ⊙ to 1.4 M ⊙ .These values encompass the range of masses and ages appropriate for our (saturation limited) ONC sample and non-cluster contaminants, based on the Besançon model (see below).The fitting routine requires a uniform grid of synthetic photometry in our filters, which we have created using Synphot (STScI Development Team 2018) and stsynphot (STScI Development Team 2020).In order to estimate the possible presence of accretion luminosity and derive mass accretion rates we use a slab model described by Manara et al. (2012) where the combination of optically thin emission generated in the preshock region and optically thick emission generated by the heated photosphere is reproduced by combining a Cloudy spectrum for a standards HII region with an 8000 K black body, respectively contributing about 1/4, and 3/4 of the total accretion luminosity (Gullbring et al. 2000).For the far-UV part of the accretion spectrum (λ ≲ 3000 Å), following France et al. (2011) we prescribe a linear decrease of flux at wavelengths shorter than the Balmer jump.
To redden the synthetic photometry, we adopt the Cardelli et al. (1989) reddening law parameterized with R V = 3.1.

Likelihood
To compare observations and simulations we use a likelihood defined as: where i runs through the different magnitudes available for each star and the labels obs and mod represent the observed and model magnitudes (mag) with their uncertainties (emag).The parameter θ instead represents the whole set of fitting parameters corresponding to the different mag i,mod .

Priors
We divide our sources in two classes of objects, each with its own specific priors: • Class A: sources having both a spectral type from Hillenbrand (1997) or Hillenbrand et al. (2013) and a distance estimate from Gaia DR3 (Gaia Collaboration et al. 2021).This class has strong priors 2 downloaded from http://svo2.cab.inta-csic.es/theory/newov2/ that allow pinpointing with high accuracy both luminosity and mass.In particular: for each star, a prior on the mass is derived from the Hillenbrand et al. (2013) temperature when available, or converting the spectral type from Hillenbrand (1997) to effective temperatures using the Luhman et al. (2003) relation.The prior is assumed to be Gaussian, centered on the mass derived from the effective temperature using the model isochrones and with standard deviation given by the mass spread corresponding to two spectroscopic sub-classes (i.e. the typical uncertainty of the spectral type reported by Hillenbrand et al. 2013).
similarly, the prior on the distance is a Gaussian centered on the reported parallax from Gaia DR3 (Gaia Collaboration et al. 2021) and standard deviation given by the corresponding parallax uncertainty • Class B: these sources, typically fainter, miss either an independent estimate of the spectral type, a Gaia parallax, or both.This is a less robust sample due to the much larger uncertainties on the temperature and/or distance.Still, also for this class, we can define priors on the basis of the results for Class A as follows: for the mass prior, if the source has no independent spectral type we start combining the overall mass distribution determined for Class A objects.Since this sample is heavily biased toward the ONC, we complement it by adding the distribution of masses obtained from the Besançon model of stellar population synthesis of the Galaxy in the direction of the ONC3 .
for the distance prior, if the distance is not provided in Gaia DR3, we take a similar approach to the one adopted for the mass.We combine the distributions of parallaxes obtained for Class A objects with the distribution of distance values obtained from the Besançon model (converted to parallax).
Given the uncertainty on the membership, ONC or Galaxy, for Class B objects we also define two other priors: one on the age and anoher on the extinction.
to create the prior on the age, we start combining the distributions obtained for Class A objects and the distribution of age values obtained from the Besançon model.
for the prior on the extinction, we take only the distributions obtained for Class A objects.We do not adopt the extinction distribution provided by the Besançon model as it does not account for the extinction caused by the OMC-1 molecular cloud, which can be both very large and non-uniform over the ONC field.
To obtain the final prior probability distribution (i.e.; the probability distribution that would link the a priori knowledge about a variable before evidence is taken into account) we applied a Gaussian kernel density estimator4 to all distributions for both Class A and B. From now on we will refer to these resulting probability density functions (PDFs) as the priors, for short.

MCMC and posterior probabilities
To derive the distributions of parameters that best describe our observations we use emcee, the Foreman-Mackey (2016) Python implementation of the Goodman & Weare (2010) affine-invariant sampler.For each given set of parameters θ, we evaluate the corresponding value of the magnitudes from the models and then we use the likelihood and prior discussed above to obtain the posterior probability distribution for θ, ensuring that the individual chains converge and are thinned to retain uncorrelated samples.Basically, this approach consists in picking different points from the sample, in every n-th step.As we are dividing these points from the overall Markov chain, the dependence becomes smaller and we can achieve a mostly independent sample.
The following convergence criteria are checked every 100 iterations: • the chain must be longer than 100 times the estimated auto-correlation time, τ f , i.e. the number of steps needed before the chain loose information about its starting properties; • the estimate of τ f has changed by less than 1% since the previous check.
As long as these criteria are met, the routine keeps running until a total of 2000 post-convergence iterations are performed, or the maximum number of steps is elapsed.
Then the routine exits the iteration loop and saves the posterior distribution of the parameters to a file before moving to the next source, till the end of the catalog.

Bayesian Analysis Products
Starting from a catalog of 3399 unique sources, we perform Bayesian SED fitting using MCMC algorithm to evaluate the star's principal parameters, i.e., mass, extinction, age, distance, and the SP acc (see Section 3.1).On average, our tests show that we are able recover the spectral type determined by Hillenbrand (1997); Hillenbrand et al. (2013) within 2 subclasses (i.e., ±200 K).
Figure 5 instead, show examples of problematic fitting, generally due to the presence of a bimodal solution in one or more parameter posterior distributions, too few filters available to properly pinpoint the right spectrum to the photometry, the presence of an excess of flux in the filters most sensitive to accretion (e.g., F336W, F656N, F658N) not reflected in the other filters, or any combination of the above.This generally produces an over or under-estimate of the spectrum for the source.
Each corner plot shows the median value, which we adopted as the best-fit estimator, as well as the 68% credible intervals of each fitted parameter's marginal posterior distribution.We define this interval so that the 68% credible interval contains 68% of the total probability, with (100 -68)/2 = 16% of the remaining probability on either side.This definition coincides with the usual 1σ interval for a Gaussian distribution.Note that the posterior distributions of our parameters are generally not Gaussian nor symmetric in most cases.The marginal distributions for the single parameters are still useful for defining the parameter credible intervals but do not capture the whole information available in the full posterior distributions, with their correlation.For detailed studies of individual sources, a complete gallery is available in the online article.

Membership selection
We will discuss in the following our approach to cluster membership selection.On the basis of our analysis, we determine ONC membership as follows.First we eliminate all the sources where the MCMC fit couldn't converge, removing 468 sources from our sample and leaving us with a total of 2931 candidates.Then, to be an ONC member: • The derived parallax of a candidate member must be compatible with 2.5 ± 0.3 mas ( 402 +57 −45 pc), i.e. 1σ from the peak of the Gaia DR3 parallax distribution for all the sources falling in our field of view.This leave us with 1278 candidate sources from the initial 2931.Taking full advantage of the posterior distribution of the parallax solution (PDPS) of each source, we evaluate the ratio of the area under the curve (AOC) of PDPS that lies within the same 2.5 ± 0.3 mas interval over the total AOC of the PDDS of each source.This ratio (that we will call area ratio or AR) provides us with an estimate of the probability that the two parallaxes are compatible.
• The derived age of a candidate member must be compatible with a recent star formation event, i.e.
≤ 10 Myr.This captures the age range of interest, as previous works (e.g., Jeffries et al. 2011;Da Rio et al. 2012;Beccari et al. 2017) show that the typical age for this cluster is ∼ 2 Myr, with an age spread at most of a few Myrs.In particular, we find that ∼ 75% of our candidates have an age ≲ 2 Myr, while if we cut the age distribution at ≲ 10 Myr, ∼ 88% of the selected sources are younger than ≲ 2 Myr and ∼ 97% are younger than ≲ 5 Myr.The age selection leads to discarding other 199 sources from our candidate sources catalog, leaving us with a candidate cluster catalog of 1079 sources.
• After visually inspecting each SED fit of these 1079 cluster members, we exclude 30 candidates due the gross inconsistency between the observed photometry and the best fit (see examples in Fig. 5), reducing the number of candidates cluster member to 1049.
This sample of 1049 sources distilled from the initial 3399 stars, i.e. ∼ 31% of the initial sample, will repre- shows a source detected in multiple filters, whose solution appears bimodal for almost all parameters and that overall fails to match the observations.The bottom row instead shows the case of a source detected in only the ACS filters that appears to be bimodal in the Av and the SPacc parameters, and that produced an overestimated spectrum of accretion (blue line) that can not be accounted for by the other filters.sent our reference catalog of candidate cluster members (or reference catalog for short).Their AR distribution shows a peak at a median value AR∼ 0.89 +0.11 −0.17 .This confirm that the majority of members in our reference catalog have high probability that their estimated distance distribution is compatible with the Gaia DR3 distance distribution for the ONC of 402 +57 −45 pc.To assess the number of stars that may have been incorrectly classified as background objects, we generated a synthetic photometry catalog of background stars using Synphot.Starting from the stellar parameters T ef f , log Z, log g and distance provided by a Besançon model of the Milky Way at the coordinates of the ONC, we determined the photometry of each source using phoenix spectra.For each filter, the magnitudes were extracted after reddening the spectrum using the sum of the extinction provided by the Besançon model and the Scandariato's ex-tinction map for the OMC (Scandariato et al. 2011), for a randomly generated position within our field of view.Comparing the number of background stars with the Besançon predictions, we are able to account for ∼ 90% of the predicted background stars in the ACS redder filters and ∼ 80% in the ACS blue and WFC-IR filters filters.Overall, this test suggest that the significant fraction of background sources determined by our approach is compatible with the predictions from the Besançon model, within the uncertainties on the extinction due to the non-uniform column density of the OMC.
As anticipated in Sec. 2, we cross-match the ACS (yellow) and WFC3 (red) surveys.Discussing the WFC3 dataset, Robberto et al. ( 2020) present a detailed Bayesian analysis to disentangle bona-fide lowmass stars and substellar cluster members from back-ground sources based on the presence of water in absorption shown by F130N-F139M color.
Figure 6 shows the overlap between ACS and WFC3 surveys and their respective characterization as candidate cluster (green) and not cluster members (blueeither background, foreground or where the sorting approach failed).It turns out that of 4504 WFC3 sources, only 2945 are also present in our ACS catalog (out of a total of 3399 ACS initial sources) or, in other words, ∼ 87% of ACS sources overlap with WFC3.We must remind the reader that in the region of overlap between ACS and WFC3 some sources can be labeled cluster in one survey and not cluster in the other, and vice-versa.Indeed, the total number of overlapping sources ( 2945) is not equal to number of overlapping candidates cluster (757) plus the number of the overlapping not candidates cluster (1237), as we need to include also the ACS candidates labeled as not cluster by WFC3 ( 225) and the WFC3 candidates labeled as not cluster by ACS (726).
From the analysis of the ACS reference catalog we find 67 sources not listed in the WFC3 catalog, while 225 sources are instead discarded by Robberto et al. (2020) on the basis of their position on the CMD.Those are recovered by our MCMC analysis thanks to a compatible solution with the Gaia DR3 estimated distance and age of the cluster.
Analyzing the WFC3 sample, out of 4504 sources we identify 1552 candidate cluster members on the basis of the Bayes Factor (BF) log(BF ) >= 0, corresponding to a ∼ 34% ratio very similar to the ∼ 31% we reported for the ACS survey.Of these 1552 sources, 278 are not listed in our initial ACS catalog.Of the remaining 1290 WFC3 candidate sources, 726 are not identified as ACS candidate sources by our method, either because the MCMC fit puts them out of the cluster (418 cases) or because it fails to converge to a solution.
Since we cannot exclude that both the 308 sources missing the MCMC fit and the 278 sources not detected by our survey (but still identified as candidate cluster by Robberto et al. 2020) are actually part of the cluster, we will add them back creating a second expanded catalog of candidate cluster members.Lacking the full distributions provided by the MCMC analysis, we will adopt the values for the mass, extinction and photometry from Robberto et al. (2013) for these extra sources, but since only 497 of the 586 missing sources have an estimate of the mass and extinction, we will add only these last sample to our expanded reference catalog, counting 1546 sources in total.
Figure 7 shows the histograms of the estimated values of each parameter for all the sources in our catalog.The mass, age, and distance panels, in particular, show a mix of low-mass (≲ 0.3 M ⊙ ) and more massive (∼ 1 M ⊙ ) stars, of very young (≲ 10 Myrs) and very old (≳ 1 Gyr) stars, as well as close ones (log(parallax) ∼ 0.4 mas corresponds to ∼ 400 pc) and far away (log(parallax) >> 0.4 mas).All these distributions are compatible with two distinct populations: one younger and closer to us, the ONC (yellow), and one older and far away (the background stars -blue).Moreover, the parallax distribution shows a small population of sources with a parallax much smaller than the one expected for ONC cluster member (red -log(parallax) > 0.4 mas) compatible with the presence of an additional population of foreground objects that are indeed expected on the line of sight.
For the candidates cluster distributions from the reference catalog, the median values of each distributions are: mass 0.21 M ⊙ ; A V 4.35 mag, age 1.07 Myr, parallax 2.52 mas, and SP acc 0.08.
Figure 8 shows the color magnitude-diagrams (CMD) for eight combinations of simultaneously observed bandpasses, i.e.ACS: m435 vs. m435-m555, m555 vs. m555-m775, m775 vs. m775-m850, and m850 vs. m775-m850, on the top row and WFPC2: m336 vs. m336-m439, m439 vs. m439-m814, m656 vs. m656-m814, NICMOS: m110 vs. m110-m130, and WFC3-IR: m130 vs. m130-m139 on the bottom row.All sources with photometry in each given pair of filters are plotted, together with our model 1 Myr and 10 Myr isochrones and the A V = 1 reddening vector, for comparison.Cluster members from the reference catalog are plotted as filled black circles, while the additional sources from the expanded catalog are marked in grey.Non-members (either because they can be tagged as field stars or due to unreliable estimates of the source's properties) are open circles.Grey circles mark instead the new sources added to create the expanded catalog as discussed above.In all these diagrams the positions of the cluster members selected using the criteria explained above are generally compatible with the locus of young stellar objects (once extinction and accretion are taken into account).The systematic departures at the bright end are due to the different saturation thresholds in the various filters.
From the masses and ages derived by the fitting routine, we have also retrieved the corresponding stellar luminosity (L ⋆ ) and effective temperature (T eff ) interpolating over our family of BT-Settl isochrones.The physical parameters estimated through our Bayesian analysis for all 3399 sources of our main catalog are provided in Table 5.The breakdown of the table is as follows: the first two columns show the ID as provided from the StraKLIP pipeline for cross-identification and the membership flag.The following columns show R.A. and Table 5. Sample of the physical parameters of all the sources included in this work.3387.8 +76.5 −87.6 −0.9 +0.3 −1.0 72.2 +0.5 Note-Table 5 is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.Note that in the published machine-readable version of this table, the errors will occupy a different column.
DEC, followed by the estimated parameters Mass, A V , Age, distance of the source from the Sun, and derived parameters T eff , logL ⋆ , R and logg with relative upper and lower error limits.

Accretion luminosity and mass accretion rate
In this Section we will discuss the process adopted to evaluate L acc and Ṁacc .
As mentioned earlier in Section 3.1, to model the intensity of the accretion spectrum in our Bayesian SED fitting we introduce an accretion parameter, SP acc , that represents the fraction of stellar bolometric luminosity that can be attributed to accretion.Having derived L ⋆ , one immediately obtains the accretion luminosity (L acc ) and therefore the mass accretion rate ( Ṁacc ) by inverting the free-fall equation that links the luminosity re-/ * >/ @ 0 * >0 @ / acc >/ @ 1 M acc >0 \U 1 @ Figure 9. Histograms of the accretion luminosity Lacc (right) and mass accretion rate Ṁacc (left) derived by the stellar parameters obtained from our Bayesian analysis (blue) with L⋆ and M⋆ overplotted on a twin axis respectively (white).... 6 is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.Note that in the published machine-readable version of this table, the errors will occupy a different column.

Note-Table
leased in when the accretion flow impacts the stellar surface with the rate of mass accretion: where R ⋆ and R in = 5 R ⋆ are the star and inner-disc radius (Gullbring et al. 1998;Hartmann et al. 1998).
In Table 6 we list the accretion properties determined for our catalog.The breakdown of this table is as follows: the first column shows the ID as provided from the StraKLIP pipeline for cross-identification.The next columns list the derived accretion parameters logL acc and log Ṁacc and relative errors for each of our sources.
Focusing on the cluster members, in Figure 9 we show the distributions of L acc and Ṁacc .In solar units, the median values are: log L acc −1.62 ± 1.43 and log Ṁacc −8.29 ± 1.39, respectively.On the same plot, using a twin axis as a reference, we also show the corresponding L ⋆ and M ⋆ distribution.
Our estimates of L acc and Ṁacc are probably affected by some selection bias inherent to the nature of the dataset we are analyzing.While it is generally true that the number of filters plays a major role in the goodness of the fit, this is particularly true for the estimate of the accretion luminosity, in particular, if the missing filters are those most sensitive to the specific characteristics of the accretion spectrum, i.e.F336W (an HST equivalent to the Johnson U band) and the F656N and F658N filters centered on the Hα line.Photometry in those bands may be missing either because the source is too faint, or even too bright with strong accretion even leading to saturation.In this last case, the other bandpasses may provide an indication of accretion, but inaccurate or underestimated values.Also, and especially in the inner regions of the ONC, the presence of circumstellar emission (e.g.proplyds) can influence the photometry in these same filters, leading to an overestimate of the accretion luminosity.Indeed, even though we removed from this analysis the known, spatially resolved proplyds, we cannot completely exclude the presence of circumstellar emission in unresolved stars.

Three-dimensional spatial distribution
The distances determined by our MCMC Bayesian routine allow reconstructing the three dimensional spatial distribution of the sources in our reference catalog.We include in the cluster sample the 30 sources with bimodal MCMC solutions discussed in Section 4.1, for which we have adopted the distances given by the highest maxima.The right panel of the same figure shows the spatial maps of the 3 populations, together with the sources with bimodal solutions.For the candidate cluster population the 2-d spatial distribution shows a remarkable concentration toward the center (identified with θ 1 Ori-C).The foreground population appear to be more uniformly spread , with a residual clustering toward the center of the cluster.One may notice that a higher concentration of foreground sources toward the nebular core is also visible in the spatial map of the Gaia foreground stars.In fact, 77% (189 sources out of 248) of our foreground sources match with Gaia foreground stars, indicating that the anomaly can be mostly driven by the Gaia prior, while 6% (16) of them are assigned to the background by Gaia, and 17% (43) are assigned to the cluster.For these 43 sources that Gaia identifies as clusters, ∼ 40% of them still have a distance compatible with the cluster within the errors provided by our fit.
The distribution of the background population shows instead a clear drop of source density along a N-S lane, a feature tracing the high-extinction caused by the main filament of the OMC1 molecular cloud (Scandariato et al. 2011).
To assess how the main parameters of the low-mass objects probed by our survey vary with the radial distance from the center, we present in Figure 11 the cumulative distribution functions (CDFs) of the number of candidate cluster members vs. distance, grouped according to their masses (left) and ages (right).The left panel in Figure 11 indicates that the more massive ( M > 0.7 M ⊙ , 185 sources, green line) have stronger tendency to cluster toward the central region than intermediate mass sources (0.2 < M < 0.7 M ⊙ , 383 sources, blue line), and these in turn are more clustered than the low-mass stars mass (M < 0.2 M ⊙ , 481 sources, gray line), that tend to be the most dispersed.In particular, 50% of the more massive sources are within about 200 arcsec (0.4 pc) vs. 290 arcsec (0.56 pc) for the intermediate mass sources 390 arcsec (0.75 pc) for the lowest mass objects.This trend indicates that mass segregation in the ONC extends to sub-solar masses, possibly down to 0.2 M ⊙ , with a degree of concentration decreasing with decreasing mass, extending into the sub-solar and possibly brown-dwarf mass range, similar to what found for masses > 1−2 M ⊙ (not probed by our survey) by Hillenbrand & Hartmann (1998).
Similarly, the right panel of Figure 11 presents the same CDF for three bins of ages.While old (> 5 Myr -green -25 sources) and intermediate age sources (between 2−5 Myr -blue -89 sources) appear similarly dispersed through the cluster, the large majority of young sources (age < 2 Myr -gray -935 sources) show a higher degree of concentration toward the central region.This result sets a rather strong upper limit to the age of the cluster.

Hertzsprung-Russell Diagram
The presence of an isochronal age spread in the HR diagram of the ONC has been previously reported and analyzed by several authors, with evidence of contamination from older populations (e.g.Prosser et al. 1994;Hillenbrand 1997;Slesnick et al. 2004;Beccari et al. 2017;Alves & Bouy 2012).In particular, Bouy et al. (2014) determined the presence of a large foreground population towards the Orion A molecular cloud, con-taining several distinct subgroups with an age spread in the range 5 − 10 Myr.Still, the fraction of foreground stars contaminating our survey area centered on the Trapezium should be small with only about ∼ 15% of the WFC3-IR sources lying in the region of significant overlap between the cluster and not-cluster members, where the Bayesian's membership analysis of Robberto et al. (2020) has larger uncertainty.The same should be true for the ACS survey given, the very similar footprint.Of the 1049 sources part of the reference catalog, only ∼ 3% of them have an age estimated between 5 to 10 Myr and ∼ 80% have a distance ≲ 400 pc, so it is quite possible that a number of them belong to a foreground population even if photometry shows that they are compatible with real cluster members.
In any case, regardless of the level of contamination by other populations, uncertainties in the measurements, presence and orientation of disk, spots, variability, anomalous extinction, etc. can all contribute to a spread in luminosity and in turn to the observed age spread (Reggiani et al. 2011).Some of these effects may be associated with photometric variability.Star spots, in particular, modulate the luminosity of a star (Gully-Santiago et al. 2017;Gangi et al. 2022) and low mass young stars show variability that can be attributed to large spots created by strong magnetic fields (Johns-Krull & Valenti 1996).Besides variability, models show that stars in the mass range 0.1 − 1.12 M ⊙ with ≳ 50% Table 7. Minimum (top) and maximum (bottom) detectable mass, in solar units, as a function of filters and extinction.Filters are separed in four blocks for the ACS, WFC3/IR, NICMOS-3 and WFPC2 cameras, in the order.
With all these caveats, we observe a peak in the age distribution with a median value at ∼ 1.1 ± 0.1 Myr.About ∼ 88% of the sources are younger than 2 Myrs, while only ∼ 3% of the total sources have an age estimate above 5 Myrs.Our results are consistent with a single major event of star formation in the ONC.
In Figure 12 we present the Hertzsprung-Russell Diagram (HRD) for the reference catalog of 1049 cluster sources in the form of Hess diagram.We bin them in boxes of ∼ 57 K by ∼ 0.076 log 10 L ⋆ and color-code their number according to the scheme presented in the figure.The 1, 3 and 10 Myr isochrones are also shown for comparison (see Section 3.1.1).
In comparison to the previous HRDs of the ONC (e.g.Da Rio et al. 2012), our new HRD appears well constrained by the 1-3 Myr isochrones, with very few sources up to 10 Myr.There is still a scatter in luminosity for effective temperatures in the range of 4000 − 3000K (i.e., the region dominated by our class A sources).When T eff decreases below ∼ 3000K, our solutions rely on fewer filters and thus the priors have an increasingly stronger weight.Besides, the range of luminosities spanned by the different isochrones becomes narrower for T eff ≲ 3000K.Both factors, the first observational and the second theoretical, contribute to the strong concentration of the sources around the 1−3 Myr isochrones.

Empirical Initial mass function
Deriving an accurate Initial Mass Function (IMF) for a rich, young stellar cluster like the ONC allows in principle to discern its variations vs. e.g., radial distance from the center, isochronal age, binarity.To achieve this goal, one must obtain a complete, extinction-limited list of reliable mass values.In the case of our dataset, the stellar masses have been estimated combining different filter sets, and within each of them one will always find two bandpasses constraining the lowest and highest mass probed, for any given exposure time, extinction, and assumed isochrone.Owing to the spread of extinction values toward the region (see Figure 7) and the non uniformity of exposure times even within a single filter due e.g. to the the overlaps between adjacent frames, it is nearly impossible to derive homogeneous values of the limiting mass and therefore an extinction limited sample from our master catalog.In practice, our data analysis strategy aimed at combining a diversity of information, Table 8.Summary of the results for the broken power law model.α0, α1 and α2 represent the low, intermediate and high mass slopes, while m lm and m mh represent the low-to-mid and mid-to-high transition masses.
This work -0.17 −0.17 Table 9. Summary of the results for the log-normal model.mc represent the characteristic mass, while σ is the width of the log-normal, m hm is the transition mass and α hm is the slope of the power-law for the high-mass regime.both priors and datasets, in order to determine the most accurate parameters of each source, is less than optimal if one wants to derive global parameters, such as the IMF, that require a homogeneous sample free from any selection bias.On the other hand, it is worth to assess our mass distribution since our methods allows retrieving sources that may be systematically rejected in more selective, controlled catalogs.To illustrate the range of masses probed by the different filter combinations, we show in the top part of Table 7 how the minimum mass reached in each filter depends on the extinction, and similarly the lower portion shows the corresponding mass at saturation limit.Building this table we have adopted our reference isochrone and the typical exposure times of the surveys, a σ ∼ 0.1 uncertainty for the faint magnitude limits, and neglected the narrow-band Hα passbands of WFPC3 and ACS as they may be strongly affected by accretion or circumstellar emission.
The table shows how high extinction values push lowmass sources increasingly out of reach, while high-mass stars become increasingly measurable as they fall below the saturation limits.Therefore, a filter set composed by our ACS bandpasses will be generally limited at the low-mass limit by the bluest filter, F435W, as it is unable to reliably detect masses M < 0.025M ⊙ even when A V = 0. Viceversa, the upper mass limit will be determined by the wideband F805LP filter, cutting off masses M > 1.345M ⊙ even with A V = 9.
Considering all these caveats on the mass, we downselected the expanded catalog of cluster members to 1428 sources, from the initial 1546 sources (see Section 4.1), and derived an ONC mass distribution as follows.Taking full advantage of our MCMC calculations, instead of adopting just the final solution for the mass of each source, we consider the last one hundred values determined by each walker as it approaches convergence.This provides us with a rather rich and robust statistical distribution of all parameters, mass in particular, compatible with the data and the priors.Then, extracting randomly for each source one value from its mass distribution, we produce a particular realization of the mass distribution in ONC.We repeat this process a hundred times, each iteration providing us with a slightly different IMF, that once averaged over the same constant bins produce the mass distribution represented by the black markers in Figure 13.The error bars indicate the 1σ spread of the data in each bin, the only anomaly being the dip observed in the bin at about 0.07 M ⊙ .This is most probably a selection effect.Stars falling in this bin are typically at the sensitivity limit of the the bluer ACS, and WFPC2 filters, yet massive enough to be poorly classified using their shallow WFC3 1.4 µm absorption band (see Section 2 for more details).The two effects combined may cause MCMC to provide degenerate or erroneous solutions.The aforementioned clustering of sources foreground objects toward the center of the cluster confirms that we do not account for all cluster mem-7 eff > .@ORJ 10 / )>/ @ 0\U 0\U 0\U 1 bers in the central region.On the other hand, our estimates -both mass and distance, i.e. membership-are more robust for larger masses, where we generally have a more complete SED, and for smaller masses where the WFC3 1.4 µm absorption band becomes prominent.The deficit of stars can be reconciled by including the stars that Da Rio et al. (2012) classified as belonging to this mass bin.This operation adds back 39 sources in the mass bin between 0.049 and 0.095 M ⊙ , obtaining the point represented by the red triangle in Figure 13 and a final catalog of 1467 total sources.
It is well established that for the low-mass stars and BDs (which is the region of the mass spectrum most relevant for our study), the IMF can be generally approximated with a shallower power law (ξ(log m) ∝ m −α ; Kroupa 2001) or a log-normal function (ξ(log m) ∝ e −(log m−log mc) 2 /2σ 2 ; Chabrier 2003).Adopting a Bayesian approach and an MCMC algorithm, we use both forms to fit our measured mass distribution in the mass range 0.025 -1.34 M ⊙ .The resulting distributions are shown in figure 14, with the best fit parameters presented in Tables 8 and 9.For reference, are also reported the values from the canonical IMFs (Kroupa -Kroupa  2012) and Gennaro & Robberto (2020).
When we compare the results for the broken powerlaw case from Table 8, our low-mass slope, α 0 = −0.17 is compatible within the 1 σ with Da Rio et al. ( 2012) measurement.The discrepancy vs. Gennaro & Robberto ( 2020) is due to our data not accounting for unresolved multiple systems with separation closer than 0.1" (or ∼ 40AU at the distance of the ONC; Strampelli et al. 2023), while Gennaro & Robberto (2020) leave the binary fraction free to vary in their bayesian model.The transition mass between the low and mid ranges, m lm = 0.20 M ⊙ , falls between the value predicted by Kroupa (2001) and the one observed by Da Rio et al. (2012), and is also in very good agreement with the prediction from Gennaro & Robberto (2020).For the intermediate-mass slope, α 1 = -1.34,we are in very good agreement with both predictions from Kroupa (2001) and Gennaro & Robberto (2020).On the other hand, we observe a shallower α 1 slope compared to the For the log-normal case (Table 9), we determine a peak mass compatible with Gennaro & Robberto (2020) IMF, but with smaller σ (our distribution is narrower).Similar to the power-law case, our derived mass function is indeed expected to be more similar to a "system" mass function as we do not try to resolve close-in companions in our data.Compared to Da Rio et al. (2012), we observe a smaller peak mass and a larger distribution (higher σ value).
Over all, for masses ≳ 0.1 M ⊙ our empirical IMF is in good agreement with a Chabrier (2003) System IMF.

Accretion rates in the ONC
In Figure 15, left panel, we present the relation between L acc and L ⋆ for our selected cluster sources.Similar to Sec. 5.2, we apply a stronger selection to reference catalog including only the sources with AR > 50%.The upper solid red line limits the area above which the majority of sources coincide with proplyds, according to the catalog by Ricci et al. (2008), making the estimate of L acc unreliable.The bottom dash-dotted red curve instead marks the stellar photospheric limit (see Eq. 1 in Manara et al. 2017a) for 1 Myr sources, below which the chromosphere is as bright or brighter than the accretion excess and therefore prevents disentangling the L acc contribution.
Selecting only the the reference catalog between these two lines, and applying the stronger cluster selection discussed above, we can build a plot showing the Ṁacc -M ⋆ relation for 631 sources (Figure 15, right panel).Similar to the cases of Chamaleon I and Lupus (Alcalá et al. 2017;Manara et al. 2017b), Ṁacc appear to increase with M ⋆ .The scatter is large, but the relation appears steeper for low-mass objects, suggesting that one may consider either a single or a broken power-law relation.We tested which approach provides a better description of our data (see Appendix A) finding that a segmented linear relation (corresponding to a broken power law in the linear scale) is preferable over a linear fit.With the already mentioned caveats about saturation and biases, we can provide for 0.01 M ⊙ ≲ M ⋆ ≲ 0. A broken power-law relation for both L acc and Ṁacc distributions is a new result compared to previous studies of accretion in young star-forming region (e.g.; Manara et al. 2017a), and in particular for this cluster (Manara et al. 2012), where a linear relation was generally preferred as there was not enough evidence to distinguish between them.For both parameters, the distributions appear strongly dependent on M ⋆ (or L ⋆ ) below a M ⋆ ≲ 0.21 M ⊙ (or L ⋆ ≲ 0.23 L ⊙ ).After, the trend becomes flatter, indicating that the Ṁacc (L acc ) remains more or less constant as M ⋆ (L ⋆ ) increase.However, there is a large spread of values (about 2 orders of magnitudes) for both distributions.
Our results on the L acc -L ⋆ can be compared with the relations obtained in other SFRs like e.g.ρ−Ophiuchi (Natta et al. 2006), σ−Orionis (Rigliaco et al. 2011), Chamaeleon I (Manara et al. 2016a, 2017b) and Lupus (Alcalá et al. 2014(Alcalá et al. , 2017)).If we only consider the region where L ⋆ ≲ 0.23 L ⊙ (before the braking point) we find a good agreement between our portion of the fit and the one obtained for Chamaeleon I, while we find a more loose match (within 3σ) with the other SFRs.Similarly, the slope of the Ṁacc − M ⋆ relationship is in broad agreement with other SFRs (e.g.Chamaleon I, Lupus and NGC 2264: Alcalá et al. 2017;Manara et al. 2017b;Venuti et al. 2014), where the data seem to better support a double power-law fit with lower mass stars having a more rapid decrease of Ṁacc compared to higher mass stars (even though the authors could not completely rule out the linear dependence).In general, when comparing results between different star forming region, age differences should be taken into account.Both the small variation observed with the ages and the high uncertainties related to age estimations itself (e.g., Soderblom et al. 2014) will make any detailed analysis very difficult to perform.
The slope of the Ṁacc relation can also be compared with theoretical models and be used as a proxy to evaluate the main processes behind the dispersal of protoplanetary disks (Clarke & Pringle 2006).Typical values suggest power laws with exponent ∼ 1.6 − 2 and spread of ∼ 1 − 2 dex.(Alcalá et al. 2014(Alcalá et al. , 2017;;Manara et al. 2016bManara et al. , 2017a;;Hartmann et al. 2016).Recent works, however, show that accretion variability alone is not enough to explain the spread observed in Ṁacc measurements (Fischer et al. 2022;Flaischlen et al. 2022) that must therefore be also related to physical processes such as disk evolution.Is still not yet clear in what measure the Ṁacc ∝ M 2 * relation is due to the disk evolutionary process or if it reflects how the initial conditions scale with stellar mass (Manara et al. 2022, and references therein).For example, Dullemond et al. (2006) were able to produce a steep Ṁacc -M ⋆ relation resulting from the imprint of the initial conditions for the formation of the star-disk system.A simple model of disk formation and evolution from collapsing cores provides a Ṁacc ∝ M ⋆ 1.8 relation, provided that cores of all mass have a similar distribution of rotation periods.
The ONC environment is characterized by strong UV flux that affects the structure and evolution of the disks closer to the central OB stars Winter & Haworth (2022).However, from our finding, at least for low-mass stars, the Ṁacc -M ⋆ relation is remarkably similar to the one observed in low-mass SFRs suggests that external photoevaporation does not play a significant role on the majority of young circumstellar disks.The slope of the relation, therefore, may be largely determined by the initial conditions set at the onset of the star formation process, which may be quite similar between regions that eventually form clusters of different sizes (Manara et al. 2022).
Photoevaporation by the disk central star can also play a role in the Ṁacc -M ⋆ slope.X-ray-driven photoevaporation produces a slope between 1.6 -1.9 (Ercolano et al. 2014) while UV-driven photoevaporation provides a value of ∼ 1.35 (Clarke & Pringle 2006).In this context, the result we obtain for M ⋆ ≲ 0.21 M ⊙ is probably more in line with the X-ray-driven scenario rather than the UV one, even though both the X-ray-driven and collapsing cores scenario provide shallower slopes compared to ours.
A consequence of a bimodal distribution on the Ṁacc -M ⋆ plane is a different evolutionary timescale for disk accretion around stars with different masses: i.e., disks around stars with M ⋆ ≲ 0.2 M ⊙ will evolve faster toward lower values of Ṁacc .The surveys in σ−Orionis (Rigliaco et al. 2011), in the same ONC (Manara et al. 2012), as well as the spectroscopic survey of L1641 (Fang et al. 2013), arrived to similar conclusions, even though we need to point out that a faster evolution for lower-mass stars is opposite to predictions by Alexander & Armitage (2006), who postulated that the viscous timescale in- A reason behind the bimodal distribution can be found in two different physical regimes at play in the early stages of disk formation and evolution (Vorobyov & Basu 2009).These authors suggested that disk selfgravity will drive large accretion torques for stars more massive than ∼ 0.2 M ⊙ soon after formation.This in turn will leave less material available in the Class II phase and therefore a lower accretion rate which will flatten the relationship.For lower mass stars instead, self-gravity is less important and disks evolve more viscously.In particular, the least-squares fit their model data provide as an exponent: n = 2.9 ± 0.5 if M ⋆ < 0.2M ⊙ n = 1.5 ± 0.1 if 0.2 ≤ M ⋆ < 3M ⊙ (5) while the corresponding fits to the observation from (Muzerolle et al. 2005, and references therein) for TTSs and BDs of age 0.5 − 3.0 Myr produce: Remarkably, the (Vorobyov & Basu 2009) model predicts a a similar relation for M ⋆ ≲ 0.21 M ⊙ compared to ours with our estimate positioning itself between the model and the observational data relations.For M ⋆ > 0.21 M ⊙ instead, our fit is much shallower, even though we have to argue here that our selected sample is limited to M ⋆ higher masses due to saturation of some or more filters (only ∼ 15% of sources in our selected catalog have higher masses).On the other hand, their data reach up to 3 M ⊙ , so we are not covering quite the same mass range.

CONCLUSION
We have revisited the Robberto et al. ( 2013) HST /ACS/WFC ONC catalog in filters F435W, F555W, F658N, F775W, and F850LP to provide updated estimates of the average magnitudes in each filter.We then combined the photometry from ACS and WFPC2 with the more recent HST /WFC3-IR measures from Robberto et al. (2020) in filter F130N and F139M.The resulting dataset contains the photometric data derived from the two HST Treasury programs dedicated to the ONC, spanning almost 15 years.We adopted a Bayesian SED fitting to derive reliable estimates of the principal stellar parameters (i.e., temperature, extinction, age, accretion, and distance) for the majority of the sources in the catalog and estimate general properties for the cluster.In particular, the three dimensional study of mass distribution for bona-fide cluster members shows that mass segregation in the ONC extends to subsolar masses, possibly down to ∼ 0.2 M ⊙ .From the age distribution we unveil a star formation history heavily peaked at ∼ 1.1 ± 0.1 Myr, with a small tail extending up to 10 Myrs.Our estimates, model dependent, do not support the idea of multiple episodes of significant star formation in the ONC.Training our dataset on the group of cluster members with known parallax, we evaluate an average distance for the cluster of 397 ± 17 pc, and similarly for the extinction distribution we estimate an median extinction of 4.35 ± 4.10 mag.
The inclusion of accretion in our fitting procedure allows us to test the relationship between L acc and Ṁacc vs. the stellar parameters.If we consider only the region where L ⋆ ≤ 0.16 L ⊙ , we also find good agreement between the L acc -L ⋆ relation for this cluster and a less massive system like e.g.Chamaeleon I, while we find a more loose agreement (within 3σ) with σ−Orionis, ρ−Ophiuchi, and Lupus SFRs.Overall, we find evidence for a bimodal relation for L acc and Ṁacc for this cluster.Both can be parameterized as broken powerlaws with slopes log L acc 1.46 ± 0.06 for log L ⋆ ≤ −0.63 and 0.42 ± 0.40 otherwise, and log Ṁacc 1.98 ± 0.10 for log M ⋆ ≤ −0.68 and 0.13 ± 0.41 otherwise.
The result we obtain for M ⋆ ≲ 0.2 M ⊙ support a scenario where the excess emission is dominated by Xray-emission from the central star rather than the UV one, while external photoevaporation doesn't look like playing a significant role on the majority of young circumstellar disks.The slope of the relation, therefore, may be largely determined by the properties of the central stars, and therefore by the initial conditions set at the onset of the star formation process, rather than by the dynamical evolution of the cluster.These conditions may be quite similar between regions that eventually form clusters of different sizes (Manara et al. 2022).
The shape of the IMF turns out to depend on the particular dataset, as the type and number of filters used to derive the source parameters is not homogeneous.This introduces selection biases that need to be carefully assessed when deriving global systems parameters such as the IMF.Overall, for masses ≳ 0.1 M ⊙ , we find our results to be compatible with a canonical Chabrier System IMF.
This work underlines the potential and limitations of multi-color photometric surveys to characterize the main physical properties of star forming regions.

Figure 1 .
Figure 1.Total area coverage of the different instruments adopted in this work.From top to bottom (left to right): ACS, WFPC2, NICMOS (superimposed to a color-composite JHK image of the Orion Nebula from 2MASS; Robberto et al. 2013) and WFC3 (superimposed to a mosaic obtained from both space and ground telescopes data; Robberto et al. 2020).For reference, the ACS mosaic covers ∼ 627 arcmin 2 , while WFC3-IR covers ∼ 486 arcmin 2 .

Figure 2 .
Figure 2. Photometric errors as a function of magnitude for the five ACS filters.

Figure 3 .
Figure3.Flowchart for the ONC sources parameter estimations.The light grey area represents the whole MCMC run while the dark one represents each single walker loop.The ensemble of all walker loops represents the simulation for one single star.The green colored blocks represent fixed ingredients for the simulation while in blue are shown the fitting parameters and in white the key steps of the Bayesian approach.The yellow block shows the stored final output of each iteration.Note that all stars are treated completely independently until the endpoint.
Fig. Set 1. ONC SED fitting Figures 4 show examples of the SED fitting (left) and corner plots (right) representing the posterior probabilities for the parameters of two random stars.The red line in the left panel shows the original spectrum evaluated for the star, while the blue line shows the slab model adopted in this work scaled by the log 10 SP acc parameter.The black line instead represents the spectrum of the star corrected by the slab model.The filters adopted for each fit are shown as circles color-coded by instrument.

Figure 4 .
Figure 4. Sample of SED fitting (left) and corner plots (right) for cluster target sources in the catalog.The red line in the left panel shows the original spectrum for the star, while the blue line shows the slab model adopted in this work scaled by the log10SPacc parameter.The black line instead represents the spectrum of the star combined with the slab model.The filters utilized for each fit are shown as circles color-coded by their respective instrument.On the right panel, the peak (blue dotted line) and the limits of the 68% credible interval (black dotted lines) are reported for each parameter.The black line in the histograms shows the Gaussian KDE.The complete figure set (1049 images) is available in the online journal.

Figure 5 .
Figure 5. Similar to Figure 4, this Figure show two examples of problematic outcome from our Bayesian fit.The top rowshows a source detected in multiple filters, whose solution appears bimodal for almost all parameters and that overall fails to match the observations.The bottom row instead shows the case of a source detected in only the ACS filters that appears to be bimodal in the Av and the SPacc parameters, and that produced an overestimated spectrum of accretion (blue line) that can not be accounted for by the other filters.

Figure 6 .Figure 7 .
Figure 6.ACS/WFC3 ONC Venn diagram showing the overlap between ACS and WFC3 surveys and their respective characterization as candidate cluster/not cluster members.

Figure 8 .
Figure 8. CMDs for ACS, WFPC2, NICMOS, and WFC3-IR filters for the whole sample of unique sources in the catalog.The dark dots mark the reference catalog sources, while the grey ones mark the the new sources added to create the expanded catalog.The red continuous and dashed-dotted line shows the 1 and 10 Myr isochrone for each CMD.Typical error bars averaged for bins of magnitude are shown on the right of each plot.

Figure 10 .
Figure 10.Spatial distribution of sources in the catalog divided by membership.Left: 3-d spatial representation of foreground (blue), candidate ONC members (black), candidate ONC members with not converging MCMC fit (gray) and background (red).An interactive 3-d version of this plot is provided in the online version of this paper.Right: 2-d view for the same groups.A yellow marker shows the position of θ 1 Ori-C on each plot.The left panel of Figure10shows the three main populations that can be discerned: the candidates background sources (red dots -1303 sources) to which the fitting procedure assigns an average distance, d ∼ 3 Kpc (log d(pc)∼ 3.6), the candidates cluster members at about 400 pc (log d(pc)∼ 2.6, black dots -1044 sources), and the population of candidates foreground sources (blue dots -350 sources) at d ∼ 250 pc (log d(pc)∼ 2.4).We include in the cluster sample the 30 sources with bimodal MCMC solutions discussed in Section 4.1, for which we have adopted the distances given by the highest maxima.The right panel of the same figure shows the spatial maps of the 3 populations, together with the sources with bimodal solutions.For the candidate cluster population the 2-d spatial distribution shows a remarkable concentration toward the center (identified with θ 1 Ori-C).The foreground population appear to be more uniformly spread , with a residual clustering toward the center of the cluster.One may notice that a higher concentration of foreground sources toward the nebular core is also visible in the spatial map of the Gaia foreground stars.In fact, 77% (189 sources out of 248) of our foreground sources match with Gaia foreground stars, indicating that the anomaly can be mostly driven

Figure 12 .
Figure 12.Hess diagram for the HRD for ONC candidate cluster sources.Luminosity and temperature are derived from the best fit of the models.Three isochrones (1, 3, 10 Myr) are shown for comparison.The colors represent the number of stars present in each box.

Figure 13 .
Figure 13.Comparison between the derived ONC mass distribution (black) and the canonical IMFs: Kroupa (yellow), Chabrier singles (blue) and systems (green).Each black dot show the average number of stars in each bin obtained form all the mass distribution histograms, with the error bars including the 1 sigma spread in each bin and the relative Poisson noise .The red marker show the correction obtained adding Da Rio et al. (2012) sources in the mass bin 0.049 -0.095 M⊙.

Figure 14 .
Figure 14.Comparison between the empirical ONC IMF (black) and from 500 draws from the posteriors (gray) and the canonical IMFs: Kroupa (left), Chabrier singles and systems (right).Also shown as reference in both plots the IMFs obtained by Da Rio et al. (2012) and Gennaro & Robberto (2020).one observed by Da Rio et al. (2012), but their large uncertainty makes the results still compatible.For the log-normal case (Table9), we determine a peak mass compatible with Gennaro & Robberto (2020) IMF, but with smaller σ (our distribution is narrower).Similar to the power-law case, our derived mass function is indeed expected to be more similar to a "system" mass function as we do not try to resolve close-in companions in our data.Compared to Da Rio et al. (2012), we observe a smaller peak mass and a larger distribution (higher σ value).Over all, for masses ≳ 0.1 M ⊙ our empirical IMF is in good agreement with a Chabrier (2003) System IMF.

Figure 15 .
Figure15.Accretion luminosity (Lacc) as a function of stellar luminosity (L⋆) for ONC sources with relative error bars.Also shown are the limits derived due to chromospheric emission for a 1 Myr old object (lower dash-dotted red curve) and the limit above which we expect the majority of sources to be proplyd (continuous red line), as explained in the text.Right: Similar to left but for the Mass accretion rate ( Ṁacc) as a function of stellar mass (M⋆) for ONC sources.The thick blue dashed segmented line represents the best fit obtained from the data using a piecewise linear function (equation A7).The colored area provides a relative visualization of the errors on the fit.creases with stellar mass to explain the Ṁacc ∝ M ⋆ 2 relation.A reason behind the bimodal distribution can be found in two different physical regimes at play in the early stages of disk formation and evolution(Vorobyov & Basu 2009).These authors suggested that disk selfgravity will drive large accretion torques for stars more massive than ∼ 0.2 M ⊙ soon after formation.This in turn will leave less material available in the Class II phase and therefore a lower accretion rate which will flatten the relationship.For lower mass stars instead, self-gravity is less important and disks evolve more viscously.In particular, the least-squares fit their model data provide as an exponent: Table 7 provides the saturation/sensitivity

Table 2 .
Number of sources detected in each filter.

Table 3 .
Sample of the photometric data available for all the sources included in this work.

Table 4 .
Range of values explored for each parameter during our MCMC analysis.From left to right: the stellar mass in Solar units, visible extinction, stellar age in Myr, the rescale factor for the slab model we use to model accretion, and the parallax (corresponding to an interval between 167 pc and 50 Kpc).

Table 6 .
Sample of the accretion properties of all the sources included in this work.