HUBBLE SPACE TELESCOPE MEASURES OF MASS ACCRETION RATES IN THE ORION NEBULA CLUSTER

As required by our method, we need to assume an isochrone that defines the intrinsic, photospheric, colors as a function of T_eff, for the ONC sources in our 2CD. As shown by Da Rio et al. (2010), current atmosphere models are still unable to lead to accurate predictions of optical colors for young (∼2 Myr) cool stars, and empirical calibrations of photospheric colors have not been specifically carried out for young PMS stars. Also, in that work the authors show that, for the same T_eff, the photospheric colors of ∼2 Myr PMS sources are not in agreement with those of MS dwarfs. Thus, we cannot either rely on calibrations of colors valid for more evolved populations. This problem therefore requires a particular care. We explain in Appendix B the procedure used to obtain the correct isochrone calibration for our 2CD, and we report in Table 1 the colors of the isochrone. It could be argued that the assumption of a single isochrone in our 2CD, given a possible age spread in Orion, may lead to false predictions. This is not the case, for two reasons: (1) the 2CD is, by its nature, independent of stellar luminosities or radii; therefore, age variations for the same T_eff can result in differences in the photospheric colors only due to changes in the stellar surface gravity; (2) in the optical, the difference in colors between PMS and MS stars can be of up to several tenths of a magnitude (Da Rio et al. 2010). This offset, scaled down to the modest age spread within the ONC (Reggiani et al. 2011), results in negligible differences in the colors versus T_eff relations.

In Figure 1, we present our calibrated isochrone, together with the observed sources. If ongoing mass accretion did not alter our observed colors, the sources would have been located on the model line, with only a spread given by extinction. Instead, a number of ONC members are displaced toward bluer colors, as expected by the presence of accretion flux excess.

For our subsequent analysis, and specifically for the derivation of bolometric stellar luminosities, we also require the bolometric corrections (BC) in the I band as a function of T_eff. These are defined as the difference between the bolometric magnitude of a star and the apparent I-band magnitude

$$\begin{eqnarray} BC &=& M_{{\rm bol}} - M_{\lambda } \nonumber \\ &=& 2.5 \log \Bigg[\frac{\int f_{\lambda } (T)\cdot S_{\lambda } \ d\lambda }{\int f_{\lambda } (T) \cdot d\lambda }\Bigg/\frac{\int f_{\lambda }(\odot)\cdot S_{\lambda } \ d\lambda }{\int f_{\lambda }(\odot) \cdot d\lambda }\Bigg]. \end{eqnarray} \tag{ 2 }$$

We computed the BCs by performing synthetic photometry on the spectra of Allard et al. (2011), and corrected the results based on the same offsets we determined in our calibration of the colors, to correct the I-band intrinsic magnitudes. The values of this empirical calibrated BC for various spectral types are also shown in Table 1. The BC_I varies between ∼0.03 and ∼ −0.70 mag in the range of temperatures of our sources (3000 K < T_eff < 8000 K).

Spectra of the flux excess produced by accretion have been modeled by Calvet & Gullbring (1998, hereafter CG98) using a hydrodynamic treatment for the infalling material and heated photosphere. It is possible, however, to reproduce a typical accretion spectrum using simpler recipes, for instance by modeling the radiative transfer within a slab of dense gas (Valenti et al. 1993) or assuming suitable combination of an optically thin emission generated in the pre-shock region and an optically thick emission generated by the heated photosphere. In the latter case, the two components can be added with a relative fraction of about 1/4 and 3/4, respectively (Gullbring et al. 2000). For simplicity, we follow this latter approach; for the optically thick component, we consider an 8000 K black body, and for the optically thin one, we modeled the emerging spectrum from a slab of gas using version 10.00 of the Cloudy software, a photoionization code last described by Ferland et al. (1998). Specifically, we assumed a slab with density n = 10⁴ cm⁻³ and temperature of ∼20,000 K. The particular choice of these two parameters has been selected after many trials, in order to obtain an accretion spectrum as similar as possible to that of CG98. For the far-UV part of the spectrum (λ ≲ 3000 Å) we follow the results of France et al. (2011) and assume a linear decrease of flux at wavelengths shortward of the Balmer continuum.

We use our photometrical colors and produce a reddening-corrected 2CD (Figure 2), assuming again the A_V values of Da Rio et al. (2012). In this way, the displacements of the ONC members from the isochrone are solely due to accretion and we can investigate this aspect, neglecting the extinction contribution. In this plot sources are represented as equally normalized two-dimensional Gaussian, with widths corresponding to the individual photometric error. This facilitates our visual inspection because, in this representation, sources with small errors appear sharp and bright while sources with poorer photometry are broader and less luminous, contributing to the plot mostly through their relative number. The sources are also color-coded according to their T_eff, highlighting that the dereddened sources are generally located at bluer colors than the photospheric ones, along stripes that depend on T_eff.

Zoom In Zoom Out Reset image size

In Figure 2, we also show the predicted displacement from the photospheric isochrone, obtained by adding a component of the basic accretion spectrum to the theoretical (calibrated) spectra. Specifically, increasing L_acc/L_tot from 0 to 1, the colors are moved along lines that start from the theoretical isochrone, at any stellar T_eff, and converge to the colors of the pure accretion spectrum, located on lower left end of the plot. For this computation, we have used the basic accretion spectrum, analogous to that of CG98. It is evident that, according to this model, accretion affects both the (B − I) and the (U − B) colors, and this prediction is somewhat compatible with the extinction-corrected photometry of some ONC members. On the other hand, there are many sources observed at (U − B) < −1, a region of this diagram which is not reachable according to this accretion spectrum model. This suggests that for these stars the Balmer jump must be higher than that of our accretion spectrum.

To correct for the discrepancy we have included another component to this spectrum in order to have a better agreement with our data. This component is modeled with an H ii region spectrum, also simulated with Cloudy, with the same density as the one previously assumed for the optically thin part but a much lower temperature (∼3000 K). This is added to the accretion spectrum in such a way that it becomes relevant only for small amounts of mass accretion. The result is shown in Figure 3: the new model predicts, consistently with our data, a larger displacement in (U − B) for low values of L_acc/L_tot. For our subsequent analysis, we will therefore consider this approach to model the color displacement due to mass accretion.

Zoom In Zoom Out Reset image size

From the modeling shown in Figure 3, we derive L_acc/L_tot for all our ONC sources. For each source, we consider its known T_eff from Da Rio et al. (2012), and therefore the unique accretion displacement curve in the 2CD that starts from isochrone point corresponding to the stellar temperature of the source. We then consider the observed photometry of the source in the 2CD and deredden it along the reddening vector until an interception with the accretion curve is found. This intersection defines, geometrically, both the amount of A_V and L_acc/L_tot for the star.

Using this method, we obtain a solution for 271 sources. In order to quantify the uncertainty in the results, we use a Monte Carlo method: for each source we displace randomly the photometry in each band within the photometric errors, assumed Gaussian. We also randomly change the stellar T_eff accounting for the uncertainty in the spectral types from spectroscopy, equal to ±1 spectral subtype (Hillenbrand 1997). We iterate this approach 1000 times for each star, deriving each time A_V and L_acc/L_tot values. From the 1000 results of A_V and L_acc/L_tot for each source, we determine the 1σ, 2σ, and 3σ intervals of the distribution of results. If an intersection is found for all the 1000 iterations, we assume our result to have a 3σ confidence; on the other hand, if for some of the iterations an intersection is not found (because, e.g., the colors lie outside the color range reachable by the accretion model plus extinction), a lower confidence level is associated with the result. We report these values in Table 2. Neglecting sources with confidence lower than 1σ, our sample comprises 245 sources. In Figure 4, we show the distribution of A_V and L_acc/L_tot for this sample. Our extinctions vary from ∼0 to ∼5, with a mean value of ∼0.98 ± 0.06 mag and only seven sources show a negative value of A_V, the lowest being A_V = −0.18  ± 0.03 mag. We assign to these sources A_V = 0 and retain them in our analysis. Comparing our results with those obtained by Da Rio et al. (2012) we obtain a very good match, measuring of a distribution of differences $\Delta A_V = A_{V_{\rm this\ work}} - A_{V_{\rm DaRio}}$ peaked at −0.03 mag with a sigma of 0.31 mag. This testifies to the accuracy of our method.

Zoom In Zoom Out Reset image size

From now on we will use these values of A_V and L_acc for all the sources for which this method gave a result. On the other hand, a small difference in the estimate of L_acc with Da Rio et al. (2010) is expected, given the higher sensitivity of our 2CD to this quantity with respect to their analysis.

The Monte Carlo approach allows to estimate the errors on our parameters. We derive individual errors from the standard deviation of the values obtained; these will be used in our subsequent analysis to estimate the errors of the $\dot{M}_{\rm acc}$.

We convert the derived L_acc/L_tot from the 2CD (Section 3.1.3) in terms of L_acc/L_☉. This is achieved considering the stellar luminosities L_* derived in Section 3.3, from the relation:

$$\begin{eqnarray} &&\frac{L_{\rm acc}}{L_{\odot }} = \frac{L_{\rm acc}/L_{\rm tot}}{1-(L_{\rm acc}/L_{\rm tot})} \cdot \frac{L_*}{L_{\odot }}. \end{eqnarray} \tag{ 4 }$$

In Table 4, we report these values for a total sample of 245 sources.

Table 4. Accretion Values of the Sources Where the 2CD Method Leads to a Determination of A_V and L_acc, According to Different Evolutionary Models

OMa	log(Lacc/L☉)	M*b	log Ageb	log $\dot{M}_{\rm acc}$b	M*c	log Agec	log $\dot{M}_{\rm acc}$c	M*d	log Aged	log $\dot{M}_{\rm acc}$d
		(M☉)	(yr)	(M☉ yr−1)	(M☉)	(yr)	(M☉ yr−1)	(M☉)	(yr)	(M☉ yr−1)
20	−3.26	0.25	7.4	−10.40	0.20	7.4	−10.30	0.23	7.4	−10.36
30	−2.95	0.17	6.9	−9.79	0.12	6.8	−9.66	0.13	6.7	−9.69
47	−2.60	0.25	6.5	−9.45	0.28	6.8	−9.49	0.25	6.5	−9.45
50	−2.42	0.29	6.4	−9.25	0.33	6.6	−9.31	0.31	6.4	−9.28
63	−3.59	0.16	6.1	−10.07	0.19	6.4	−10.15	0.14	5.9	−10.00
69	−3.17	0.17	7.0	−10.06	0.13	6.9	−9.94	0.14	6.9	−9.98
70	−4.13	0.18	7.0	−11.04	0.14	7.0	−10.94	0.14	6.8	−10.92
77	−1.45	0.40	6.3	−8.33	0.42	6.5	−8.36	0.43	6.4	−8.36
87	−2.99	0.17	7.1	−9.91	0.13	7.0	−9.77	0.14	6.9	−9.83
146	−2.31	0.32	5.7	−8.92	0.47	6.1	−9.08	0.47	6.0	−9.09
167	−1.56	0.27	6.4	−8.37	0.31	6.6	−8.42	0.28	6.4	−8.38
...	...	...	...	...	...	...	...	...	...	...

Notes. ^aOrion Master Catalog entry number (Robberto et al. 2012). ^bObtained with D'Antona & Mazzitelli (1994) models. ^cObtained with Siess et al. (2000) models. ^dObtained with Palla & Stahler (1999) models.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

There are several empirical relations between the L_Hα and the total accretion luminosity. Since, however, we have an independent estimate of the L_acc for some of our ONC sources from the 2CD, we rederive this transformation based on our data. To this purpose, we isolate the sources where both L_acc obtained with the 2CD and L_Hα are available. In particular, limiting to sources with at least a 2σ confidence (see Section 3.1.3), we isolate a sample of 148 stars. Figure 5 shows the relation between these two quantities; to account for the overall uncertainties, including the intrinsic scatter in the values in the Hα flux due to variability (e.g., Murphy et al. 2011), we add quadratically to all the estimated errors on log (L_acc/L_☉) a minimum value of 0.2, and use these combined errors for the linear regression. Our fit provides the following relation:

$$\begin{eqnarray} &&\log \left(\frac{L_{{\rm acc}}}{L_\odot } \right) = (1.31\pm 0.03) \cdot \log \left(\frac{L_{{\rm H}\alpha }}{L_\odot } \right) + (2.63\pm 0.13). \nonumber \\ \end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Herczeg & Hillenbrand (2008), fitting the Balmer Jump and veiling features and correlating those estimated for L_acc to the spectroscopical estimated EW_Hα, found a relation with similar slope (1.20 ± 0.16). Dahm (2008) measured a slope of (1.18 ± 0.26) using 14 members of the Taurus–Auriga complex, combining the previous L_acc estimate with their spectroscopical EW_Hα. De Marchi et al. (2010) assumed a unitary slope, implying the proportionality between L_Hα and L_acc, and derived another relation using the same data of Dahm. We rely on our estimate, and we thus use the relation in Equation (5) to obtain L_acc for 528 additional sources with no estimate of this quantity from the 2CD.

The completeness of these samples depends on two factors: the availability of T_eff from Da Rio et al. (2012) and the detection of the sources in U or Hα (within the WFPC2 field of view). Whereas the U-band photometry tends to be shallow for highly reddened stars, the Hα photometry is much less affected by extinction, and extends well into the BD mass range. Therefore, source detection is not expected to alter the representativeness of our stellar sample. On the other hand, the available T_eff spectral types from Da Rio et al. (2012) are characterized by a uniform and close to 100% completeness down to the H-burning limit. Therefore, we expect our sample of $\dot{M}_{\rm acc}$ to be representative of the ONC stellar population down to ∼0.1 M_☉, and the lower number of sources in our catalog compared to Da Rio et al. (2012) is mainly due to our smaller field of view. To test this hypothesis, and rule out significant selection effects in our sample of $\dot{M}_{\rm acc}$ estimates, we consider the full catalog of available stellar parameters from Da Rio et al. (2012), limited to the same (smaller) field of view of our WFPC2 survey. Then we compute, as a function of mass and age, the fraction of these sources included in our $\dot{M}_{\rm acc}$ sample. The result, obtained assuming D'Antona & Mazzitelli models, is shown in Figure 8. The error bars represent the Poissonian uncertainty from the stellar numbers in both samples. As supposed, we do not detect significant trends in the fraction of accreting sources with respect to the stellar paramters, except for a significant lack of sources below 0.1 M_☉. Above this mass, about two-thirds of the ONC members have an $\dot{M}_{\rm acc}$ measurement, regardless their mass or age.

Zoom In Zoom Out Reset image size

To investigate the variations of $\dot{M}_{\rm acc}$ with respect to the stellar age and mass, we consider all of these three quantities in a single 3D space. To this purpose, we neglect sources with suspicious ages, specifically, too young (log  (t/yr) < 5.5) or too old (log  (t/yr) >7.3).

We perform a minimum χ²_red fit, accounting for errors on the single values, assuming two different forms: (1) a simple plane (log $\dot{M}_{\rm acc}$ = A  +  B  ·  log t  +  C  ·  log M_*) and (2) a warped surface (log $\dot{M}_{\rm acc}$= A  +  B  ·  log t  +  C  ·  log M_*  +  D  ·  log t  ·  log M_*). Specifically, the first form is the simplest function to describe correlations of $\dot{M}_{\rm acc}$ with the stellar parameters; the second allows for a mixed term that accounts for a possible different evolution timescale of $\dot{M}_{\rm acc}$ with stellar mass. To perform the fit we use a standard Levenberg–Marquard nonlinear regression. In the two cases the resulting χ²_red are comparable (∼7–9), and the second one is slightly smaller. This decrease of χ²_red is not surprising, given that the first function is a subcase of the second, the latter having an additional parameter. Therefore, in order to quantitatively decide whether the most complex model is truly more representative of the real distribution of the data, or whether the simplest model is good enough given the uncertainties, we run a statistical F test. This test compares the relative increase in the residual sum of squares obtained by reducing the complexity of the model, with the relative increase of the degrees of freedom. Then, the test derives the probability that the simpler model is suitable enough. In any case, we find this probability to be lower than 0.1%, meaning that introducing the fourth parameter D in our model (i.e., a warped surface instead of a flat plane) leads to a more representative match with our data, and this is not a statistical casualness. Our fitted surface is shown in Figure 11 and described by the following equation:

$$\begin{eqnarray} &&\log \dot{M}_{\rm acc}= A + B \cdot \log t + C \cdot \log M_* + D \cdot \log t \cdot \log M_*, \nonumber \\ \end{eqnarray} \tag{ 6 }$$

where the parameters' values are reported in Table 6, also in the cases where we consider only sources with results obtained from the 2CD or with accretion estimated using the Hα luminosity.

Zoom In Zoom Out Reset image size

These relations imply that both parameters together play a role in the evolution of mass accretion. To better understand this behavior, and in particular the role of the "mixed" term described by the parameter D, we separately investigate the dependence of $\dot{M}_{\rm acc}$ on age or mass; also, in every case, we divide the sample in different ranges for these parameters. This is illustrated in Table 7 for the dependence on age, and in Table 8 for the dependence on mass, and the plots are shown in Figure 12.

Zoom In Zoom Out Reset image size

If we consider only the dependence of $\dot{M}_{\rm acc}$ on stellar ages ($\dot{M}_{\rm acc}$ ∝t^−η), our fit implies that a similar slope as that proposed by Hartmann et al. (1998) (η = 1.5) is obtained for M_* ∼ 0.13 M_☉. For higher mass stars we find decreasingly lower values of η with mass. Table 7, in particular, shows that the values of η decrease monotonically regardless of considering all sources, only those with results obtained from the 2CD or with the Hα luminosity. Thus, our data indicate that there are different accretion evolutionary timescales for sources with different mass. Surprisingly, for sources with M_* ≳ 0.5 M_☉ the slope is smaller than 1. This value is not compatible with the similarity solution model framework (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). Indeed, if we suppose that the disk evolves with time following a simple power law (M_disk∝t^−α), we have that $\dot{M} = dM_d/dt \propto t^{-(\alpha +1)}$; thus η = α + 1 > 1 in order to have a disk mass decreasing with time.

On the other hand, according to the self-similar models, $\dot{M}_{\rm acc}$ scales with time as a simple power law only in the asymptotic regime, when the age of the system t ≫ t_ν, where t_ν is the viscous timescale of the disk (set by the initial conditions and by the relevant viscosity). At earlier times, the relationship between $\dot{M}_{\rm acc}$ and t is indeed shallower. Our results might thus be explained within the framework of the self-similar solutions if higher mass stars are not yet in the asymptotic regime, implying that the viscous time is a growing factor of the stellar mass. Note that this trend is the opposite of what suggested by Alexander & Armitage (2006) in the BD regime to explain the M_* − $\dot{M}_{\rm acc}$ correlation.

Alternatively, one might question the validity of the self-similar models themselves in describing, from a statistical point of view, the evolution of protostellar disks, and in particular the assumption that viscosity scales with radius as a simple power law. Indeed, different physical processes are expected to redistribute angular momentum within the disk at different radii, with the magneto-rotational instability (MRI; Balbus & Hawley 1992) dominating at small distances from the star, the gravitational instability (Lodato & Rice 2004) dominating in the colder outer parts, and with the possibility of having extended dead zones (Gammie 1996) in between the two.

For what concerns the dependence of $\dot{M}_{\rm acc}$ over mass (see Table 8) modeled as $\dot{M}_{\rm acc}$ ∝ M^b_*, the trend is similar regardless on the method used to estimate $\dot{M}_{\rm acc}$. In particular, the slope in higher, thus the relation is steeper, when considering older ages. This also shows that at older ages more massive sources maintain a higher accretion rate.

Figure 12 illustrates once more these trends. We plot the evolution of our $\dot{M}_{\rm acc}$ values against age and M_* for different ranges of the parameters. In this figure, the slopes of the overplotted lines are obtained from Equation (6), substituting to the projected parameter (mass for the left panels, age for the right ones), the mean value within the subsample used for each particular panel.

The slower temporal decay of $\dot{M}_{\rm acc}$ that we measure in the intermediate-mass range (M_* > 0.4 M_☉) is in agreement with estimates of mass accretion rates from Hα photometry obtained in the Magellanic Clouds (De Marchi et al. 2010, 2011; Spezzi et al. 2012) and in the massive galactic cluster NGC 3603 (Beccari et al. 2010). Indeed, in all these studies, the investigated stellar samples do not reach the low-mass regime, due to poor photometric sensitivity. This suggests, as inferred by the authors of these works, that their relatively high values of $\dot{M}_{\rm acc}$ measured for isochronal ages of ∼10 Myr compared to galactic low-mass stars indicate a selection effect because of higher masses.

Our derivation of $\dot{M}_{\rm acc}$ using the 2CD may be, in principle, affected by some modest selection effects. This, as mentioned, is located in the fact that for high T_eff sources, the color displacement due to accretion is very small, limiting the possibility to identify weak accretors among intermediate mass stars. Specifically, Figure 3 shows that for (B − I) ≲ 2, the 1% L_acc/L_tot line is undistinguishable from the zero accretion one. This is simply because stars with high T_eff emit a non-negligible fraction of their flux in the U band, "hiding" a possible small accretion excess. As mentioned, however, this problem is highly mitigated since we have also derived L_acc from the L_Hα excess for all the sources for which $\dot{M}_{\rm acc}$ is not derived from the 2CD. In the end, with the Hα measurements, and even considering upper limits (see, e.g., Figure 10) the relations we find do not appear to suffer from this mentioned selection effect.

As already introduced in Section 4.2, the accuracy of the isochronal ages to represent the true stellar ages has been recently questioned. Baraffe et al. (2009) and Baraffe & Chabrier (2010) claimed that the observed spread of L_* in the HRD of star forming regions can be explained as a consequence of the protostellar phase accretion history. In particular, episodic events of intense accretion ($\dot{M}_{\rm acc}$   ∼  10⁻⁴ M_☉ yr⁻¹) during this phase can produce a large spread of values of stellar radii for stars of identical age and mass. Their results are based on the assumption of "cold" accretion and initial masses M_i between 1 M_J and 0.1 M_☉. On the other hand, Hosokawa et al. (2011) demonstrate that the apparent spread produced by different accretion histories may be important only for sources with T_eff ≳ 3500 K (∼30% of our sources). In their analysis they assume a fixed M_i = 0.01 M_☉, and this could suggest that the results obtained should be revised changing the value of this parameter. Thus, from a theoretical point of view the debate is still ongoing and a unified picture is not yet reached. At the present moment, therefore, we can only use the age inferred from the evolutionary tracks used in the past (e.g., D'Antona & Mazzitelli 1994), waiting for new isochrones based on episodic accretion models.

We stress that if the isochronal ages we assume are inaccurate, because of the aforementioned effects, this results in an overall "horizontal" broadening of our results in the $\dot{M}_{\rm acc}$ versus age plane, but the overall trends we measure cannot be strongly influenced.