PHOTOMETRIC DETERMINATION OF THE MASS ACCRETION RATES OF PRE-MAIN-SEQUENCE STARS. I. METHOD AND APPLICATION TO THE SN 1987A FIELD^*

In the current star formation paradigm, low-mass stars grow in mass over time through accretion of matter from a circumstellar disk (e.g., Lynden-Bell & Pringle 1974; Appenzeller & Mundt 1989; Bertout 1989). This accretion process is believed to occur via the magnetic field lines that connect the star to its disk and that act as funnels for the gas (e.g., Königl 1991; Shu et al. 1994). The strong excess emission observed in most of these pre-main-sequence (PMS) stars is believed to originate through gravitational energy, released by infalling matter, that ionizes and excites the gas. Thus, the excess luminosity of these objects (the classical T-Tauri stars) can be used to measure the mass accretion rate.

A reliable measurement of the rate of mass accretion onto PMS stars is of paramount importance for understanding the evolution of both the stars and their disks (see, e.g., review by Calvet et al. 2000). Of particular interest is determining how the mass accretion rate changes with time as a star approaches its main sequence (MS), whether it depends on the final mass of the forming star, and whether it is affected by the chemical composition and density of its parent molecular cloud or by the proximity of hot, massive stars.

Ground-based spectroscopic studies of nearby young star-forming regions (e.g., in Taurus, Auriga, Ophiuchus, Orion) show that the mass accretion rate appears to decrease steadily with time, from ∼10⁻⁸ M_☉ yr⁻¹ at ages of ∼1 Myr to <10⁻⁹ M_☉ yr⁻¹ at ∼10 Myr (Muzerolle et al. 2000; Sicilia-Aguilar et al. 2005, 2006). While at face value this is in line with the expected evolution of viscous disks (Hartmann et al. 1998), the scatter on the data is very large, exceeding 2 dex at any given age. This is at least in part explained by the recent finding that the mass accretion rate seems also to depend quite steeply on the mass of the forming star. Muzerolle et al. (2003, 2005), Natta et al. (2004, 2006), White & Hillenbrand (2004), and Sicilia-Aguilar et al. (2006) consistently report that the accretion rate $\dot{M}_{\rm acc}$ decreases with the stellar mass M as $\dot{M}_{\rm acc} \propto M^2$, albeit with large uncertainties on the actual value of the power-law index. On the other hand, Clarke & Pringle (2006) warn that such a steep decline might be spurious and strongly driven by detection/selection thresholds.

While observational uncertainties and potential systematic errors may well contribute to the present uncertainty, the true limitation in this field of research currently comes from the paucity of available measurements. Indeed, all the results so far obtained are based on the mass accretion rates of a small number of stars (just over 100), all located in nearby Galactic star-forming regions, covering a very limited range of ages and no appreciable range of metallicity (all clouds mentioned above having essentially solar metallicity; e.g., Padgett 1996).

The origin of this limitation is to be found in the technique used to measure PMS mass accretion rates. Regardless as to whether the latter are derived from the analysis of veiling in the photosperic absorption lines (usually at ultraviolet and optical wavelengths) or through a detailed study of the profile and intensity of emission lines (usually Hα, Paβ, or Brγ), these techniques typically require medium- to high-resolution spectroscopy for each individual object. Even with modern multi-object spectrographs at the largest ground-based telescopes, this approach is not very efficient and cannot reach beyond nearby star-forming regions or the Orion cluster, as crowding and sensitivity issues rule out but the brightest objects in more distant regions like the center of the Milky Way or the Magellanic Clouds. While in principle the PMS mass accretion rate can be estimated from the U-band excess if the spectral type is known (see, e.g., Gullbring et al. 1998; Romaniello et al. 2004), this method is potentially subject to large uncertainties in the extinction and in the determination of the correct spectral type and effective temperature when spectroscopy is not available.

In order to overcome these limitations and significantly extend the sample of stars with a measured mass accretion rate, we have developed and successfully tested a new method to reliably measure $\dot{M}_{\rm acc}$ that does not require spectroscopy. The method, described in the following sections, combines broadband V and I photometry with narrowband Hα imaging and allows us to identify all stars with excess Hα emission and to derive from it the accretion luminosity L_acc and hence $\dot{M}_{\rm acc}$ for hundreds of objects simultaneously.

The paper is organized as follows. In Section 2, we address the identification of PMS stars via their color excess and in Section 3 we show how to reliably convert the measured excess into Hα luminosity and equivalent width. Section 4 shows how the accretion luminosity and mass accretion rate can be obtained from the Hα luminosity and Section 5 presents an application of the method to the field of SN 1987A, in the Large Magellanic Cloud (LMC). A thorough discussion of all systematic uncertainties involved in the method is provided in Section 6, while a general discussion and conclusions follow in Section 7. The Appendix provides additional figures and tabular material.

One of the characteristic signatures of the accretion process on to PMS stars is the presence of excess emission in Hα (see Calvet et al. 2000). Although surveys to search for stars with Hα excess emission over large sky areas have been traditionally carried out via slitless spectroscopy (e.g., Herbig 1957; Kohoutek & Wehmeyer 1999), extensive studies in which narrowband Hα imaging is combined with broadband photometry also exist (e.g., Parker et al. 2005; Drew et al. 2005).

Recently, Romaniello (1998), Panagia et al. (2000), and Romaniello et al. (2006) used HST/WFPC2 photometry in Hα (F656N) and the R band (F675W) to provide the first direct detection of extragalactic PMS stars with strong Hα emission in the LMC. Their approach is based on the widely adopted use of the R-band magnitude as an indicator of the level of the photospheric continuum near the Hα line. Therefore, stars with strong Hα emission will have a large R − Hα color. From the analysis of the specific F656N and F675W filter response curves of the WFPC2 instrument, Panagia et al. (2000) conclude that an equivalent width (W_eq) of the Hα emission line in excess of 8 Å will result in a color excess R − Hα ≳ 0.3.

This way of identifying PMS stars is more accurate and reliable than the simple classification based on the position of the objects in the Hertzsprung–Russell (H–R) diagram, i.e., stars placed well above the MS (e.g., Gilmozzi et al. 1994; Hunter et al. 1995; Nota et al. 2006; Gouliermis et al. 2007). The reason is twofold: (1) contamination by older field stars and the effects of differential extinction may lead to an overestimate of the actual number of candidate young, low-mass PMS stars, i.e., those "above" the MS and (2) older PMS stars, already close to their MS, cannot be properly identified with this method, thereby resulting in an underestimate of the number of stars that are still forming. When PMS stars are identified by means of their R − Hα color excess, both problems disappear.

Similarly, this approach is more practical, efficient, and economical than slitless spectroscopy for studying PMS stars in dense star-forming regions in the Milky Way and nearby galaxies, which are today accessible thanks to high-resolution observations, in particular those made with the Hubble Space Telescope (HST).

We will show in Section 3 that the level of the photospheric continuum near Hα can be directly derived from the observed V and I magnitudes. This means that PMS stars with an Hα excess can be identified without actually measuring the canonical R − Hα index. On the other hand, since R − Hα has traditionally been used in PMS stars studies, we first show here how one can derive analytically the R magnitude of an object from the measured V and I magnitudes. In this way, we intend to make it possible for those who do not dispose of R-band photometry for their fields to derive it in an accurate way by using measurements in the neighboring V and I bands. The latter are more commonly used than the R band in photometric studies of resolved stellar populations and are particularly popular among HST observers.

2.1. Deriving R from V and I

The absence of remarkable spectral features between 5000 Å and 9000 Å in the stellar continuum of dwarfs and giants with effective temperatures in the range 4000 K ≲T_eff≲ 10,000 K implies that their V − R and V − I colors depend primarily on the effective temperature of the stars and less on their metallicity. In effect, V − I remains a useful color index for temperature determinations (von Braun et al. 1998; Bessell et al. 1998).

The color–color sequences in the UBVRIJKLMN passbands published by Johnson (1966) revealed a very tight correlation between the mean V − R and V − I colors of local dwarfs and giants (see Figure 1). The correlation is remarkably linear in the range 0 < V − I < 2 (i.e., 4000 K ≲T_eff≲ 10,000 K), where (V − R) = 0.58 × (V − I) (dashed line in Figure 1), and it remains linear for 2 < V − I < 4 (i.e., 3000 K ≲T_eff≲ 4000 K), albeit with a shallower slope (0.45; dot-dashed line). It is, therefore, possible and meaningful to attempt to derive the value of the R magnitude from the knowledge of V and I, particularly because no extrapolation beyond the available wavelength range is involved. The simple inversion of the relationship above gives R = 0.42 V + 0.58 I for the system of Johnson.⁶

Zoom In Zoom Out Reset image size

Interestingly, this relationship is very similar to the one that we would obtain if we assumed that the proportion with which the V and I magnitudes contribute to the R magnitude is simply dictated by the ratio of the effective wavelengths of the filters. The effective wavelengths of the V, R, and I bands in Johnson's system are, respectively, λ^V_e = 5500 Å, $\lambda _{\rm e}^{\rm }R =7000$ Å, and λ^V_e = 9000 Å, which imply that the V band should contribute for (7000 − 5500)/(9000 − 5500) ≃ 0.43 and the I band for (9000 − 7000)/(9000 − 5500) ≃ 0.57 to the R-band magnitude. These values are remarkably close to those observed.

The apparent disarming simplicity of this conclusion, however, should not lead one to think that this is the case for any photometric system. In particular, for spectral bands much wider than Johnson's, such as some of the HST/WFPC2 and HST/Advanced Camera for Surveys (ACS) filters, the matter becomes more complicated because the true effective wavelength can depend significantly on the spectral properties of the source (i.e., its temperature). In these cases, the relationship should be derived via direct observations in the specific bands or through synthetic photometry, when the properties of the photometric system are well characterized.

Due to the large number of inter-related observing modes offered by the scientific instruments on-board the HST, synthetic photometry has long been established as the most practical and efficient approach to their calibration (see Koornneef et al. 1986; Horne 1988; Sirianni et al. 2005). For this reason, reference tables exist that accurately describe the main optical components of the telescope and instruments as well as the response and quantum efficiency of the detectors as a function of wavelength. Therefore, thanks also to the availability of standardized software packages (e.g., Synphot; see Laidler et al. 2008), it is relatively easy to derive synthetic magnitudes from observed or theoretical stellar spectra for various HST instrumental configurations.

We have used both observed and theoretical stellar spectra, taken respectively from the Bruzual–Persson–Gunn–Stryker (BPGS) Spectrophotometry Atlas (see Gunn & Stryker 1983) and from the stellar atmosphere models of Bessell et al. (1998), to compute color relationships between R-like magnitudes and V − I-like colors for various HST instruments and ground photometric systems. Details on how the color transformations were derived and a table of coefficients (Table 3) can be found in the Appendix, while a brief comparison with HST data is given in the following section.

2.2. Comparison with the Data

To verify the validity of the predicted color relationships, we have used WFPC2 photometry of a field around SN 1987A studied by Romaniello (1998), Panagia et al. (2000), and Romaniello et al. (2002). The observations cover a field of 9.16 arcmin² and were collected over three epochs, namely, 1994 September, 1996 February, and 1997 July. The advantage of this catalog is that the magnitude of each star in it is individually corrected for the effects of interstellar extinction by using the information provided simultaneously by all broadband filters (see Romaniello et al. 2002 for more details). While it would be possible to include the effects of reddening in the calculations of the color relationships of Table 3, it is more appropriate in this case to correct the observed magnitudes, since considerable differential extinction is likely to be present in this field (see Panagia et al. 2000).

From the catalog of Romaniello et al. (2002), we have selected all those stars whose mean error δ₃ in the three bands F555W, F675W, and F814W does not exceed 0.05 mag, where

$$\begin{equation} \delta _3 = \sqrt{\frac{\delta ^2_{555}+\delta ^2_{675}+\delta ^2_{814}}{3}} \end{equation} \tag{ 1 }$$

and δ₅₅₅, δ₆₇₅, and δ₈₁₄ are the photometric uncertainties in each individual band after correction for interstellar extinction. A total of 9635 stars satisfy the condition set by Equation (1), out of 21,955 objects in the complete catalog.

The color–color diagram for the stars selected in this way is shown in Figure 2. The dashed line in that figure corresponds to the color relationship derived for the bands in question from the models of Bessell et al. (1998; see the Appendix) and is practically indistinguishable from the best linear regression fit to the data, which has a slope of 0.645 ± 0.003. Although these models were selected to match the metallicity of the LMC, having used the relationship derived with the BPGS atlas or that for [M/H] = −1 would have still resulted in a very good fit, within the observational uncertainties. This result confirms that, at least over the range explored here, metallicity does not play a dominant role in the color relationships.

Zoom In Zoom Out Reset image size

We provide a further proof of the reliability and accuracy of the interpolated R-band magnitudes in Figure 3, where we compare the identification of PMS stars based on the measured and interpolated R − Hα color. To distinguish the observed from interpolated R-band magnitude, we indicate the latter with R*, as explained below.⁷ In this case, from the photometry of Romaniello et al. (2002) already corrected for reddening, we have further restricted the selection to objects with δ₃ ⩽ 0.02, comprising 3241 stars.

Zoom In Zoom Out Reset image size

Diamonds in Figure 3 represent the average R − Hα color of the population and, as such, define the reference with respect to which one should look for excess emission. The dashed lines represent the theoretical color relationship obtained using the Bessell et al. (1998) model atmospheres mentioned above for normal MS stars and are in rather good agreement with the average observed colors in the range −0.1 ≲ V − I ≲ 1 (filled diamonds). The rms deviation between the model and the data over this range amounts to 0.03 mag and is dominated by the systematic departure in the range 0.2 ≲ V − I ≲ 0.7. The latter could be due to a number of reasons, including the inability of the models to accurately describe the properties of the Hα line since they are meant to be used for broadband photometry. The departure of the models from the data at V − I>1 stems from both small number statistics and the fact that the bulk of the population there comes from red giants, which have lower surface gravity and possibly lower metallicity than those assumed in the models. In the following, we will concentrate our analysis to the range −0.1 ≲ V − I ≲ 1 marked by the filled diamonds.

Figure 3(a) is based on the observed data and we indicate as large dots all stars whose observed R − Hα color exceeds the local average by at least four times the individual photometric uncertainty in the color, i.e., the one specific to each star. Consequently, these objects must be regarded as having a bona fide Hα excess above the 4σ level. There are 199 such objects in panel (a). In panel (b), we have replaced the observed R magnitude with its interpolated value R*. All but 16 of the 199 stars with 4σ excess in panel (a) have a 4σ excess also in panel (b) and are indicated as large dots. The remaining 16 objects, marked as triangles, fail this requirement as their color is closer to the local average color than four times their specific combined photometric uncertainty.

In panels (c) and (d), we study the uncertainty in the other direction, i.e., from the interpolated to the observed plane. Panel (c) simulates the case in which, no R-band data being available, one uses the interpolated R* value from Table 3 to identify stars with Hα excess emission. Stars indicated as large dots in panel (c) are those with interpolated R* − Hα color exceeding the local average by at least four times the individual photometric uncertainty in the color, i.e., those with Hα excess above the 4σ level. A total of 194 stars meet this condition. All but 11 of them would have met the 4σ excess condition if we had used the actually observed R-band magnitude. They are indicated as large dots in panel (d), while the 11 objects failing the condition are marked as triangles.

From Figure 3, we can draw some quantitative conclusions on the statistical significance of the detection of the color excess in Hα. First of all, it is not surprising that some objects detected as bona fide Hα excess stars in panel (a) are rejected in panel (b), since the typical photometric uncertainty on R* is a factor of $\sqrt{2}$ higher than that on R, as one must combine the uncertainties on V and I that concur to the estimate of R* (although the major contribution to the uncertainty on the R − Hα color of a star comes from the Hα band). For the same reason, it is to be expected that panel (c) will yield less objects with bona fide Hα excess than panel (a).

Considering the number of detections and mismatches of Figure 3, we estimate that using the interpolated R* in place of R a fraction of (199 − 194)/199 = 2.5% of the objects with Hα excess at the 4σ level would be missing and that for about 11/194 = 5.7% of the detections the assigned statistical significance would be higher than indicated by direct R measurements. We conclude that our interpolation scheme is quite reliable in that it can provide the correct identification of PMS stars in no less 94% of the cases.

While helpful to accurately identify PMS stars, the R − Hα color excess alone does not immediately provide an absolute measure of the Hα luminosity L(Hα) nor of the equivalent width W_eq(Hα), which are necessary to arrive at the mass accretion rate. In order to measure L(Hα), one needs a solid estimate of the stellar spectrum in the Hα band (i.e., without the contribution of the emission). Similarly, to measure W_eq(Hα), a solid estimate of the level of the stellar continuum inside the Hα band is needed. Since the R band is over an order of magnitude wider than the Hα filter, it cannot provide this information accurately, but only gives a rough estimate of the continuum level. In this section, we present a simple method that allows us to measure L(Hα) for the objects with an Hα excess and to derive the level of their spectral continuum inside the Hα band, determining in this way also W_eq(Hα).

3.1. The Hα Line Luminosity

The method is based on the simple consideration that, at any time, the largest majority of stars of a given effective temperature T_eff in a stellar population will have no excess Hα emission. Therefore, for stars of that effective temperature, the average value of a color index involving Hα (for instance V − Hα, and thus not just R − Hα) effectively defines a spectral reference template with respect to which the Hα color excess should be sought. Note that this is not only true for populations comprising both young and old stars, such as those typical of extragalactic stellar fields, but also for very young populations since PMS objects show large variations in their Hα emission over hours or days (e.g., Fernandez et al. 1995; Smith et al. 1999; Alencar et al. 2001), with only about one-third of them at any given time being active Hα emitters above a W_eq(Hα) ≃ 10 Å threshold (Panagia et al. 2000).

To better clarify the working of our method, we display in Figure 4 the V − Hα color as a function of V − I for a set of stars taken from the catalog of Romaniello et al. (2002). This time, we have selected a total of 4156 stars, namely, all objects with δ₃ ⩽ 0.1 mag, where the mean uncertainty δ₃ is defined as in Equation (1), but for the V, I, and Hα bands instead of V, I, and R. Therefore, no R-band information is used in this case and, in fact, no R-band data are needed. It should be noted that δ₃ is dominated by the uncertainty on the Hα magnitude, while the median value of the uncertainty in the other two bands is δ_V = 0.014 and δ_I = 0.016, respectively.

Zoom In Zoom Out Reset image size

The thick dashed line in Figure 4 represents the average V − Hα color obtained as the running median with a box-car size of 100 points. The thin solid line shows the colors in these filters for the model atmospheres of Bessell et al. (1998) mentioned above (note again the <0.1 mag discrepancy between models and observations around V − I ≃ 0.6, most likely due to the coarse spectral sampling of the Hα line in the models). At the extremes of the distribution (V − I < −0.1 and V − I>1) the density of observed points decreases considerably. To determine the reference template there, we have reduced the box-car size to 10 points and averaged the running median obtained in this way with the theoretical models (thin solid line). The thick dot-dashed line in Figure 4 shows the best linear fit to the resulting average (which is extrapolated for V − I>2).

A total of 189 objects, indicated by large dots in Figure 4, have a V − Hα index exceeding that of the reference template at the same V − I color by more than four times the uncertainty on their V − Hα values. These are the objects with a V − Hα excess at the 4σ level and, as we shall see later, most of them are bona fide PMS stars.

Since the contribution of the Hα line to the V magnitude is completely negligible, the magnitude ΔHα corresponding to the excess emission is simply:

$$\begin{equation} \Delta {\rm H}\alpha = (V-{\rm H}\alpha)^{\rm obs} - (V-{\rm H}\alpha)^{\rm ref}, \end{equation} \tag{ 2 }$$

where the superscript "obs" refers to the observations and "ref" to the reference template. Once ΔHα is determined in this way, the Hα emission line luminosity L(Hα) can be immediately obtained from the photometric zero point (ZP) and absolute sensitivity of the instrumental setup and from the distance to the sources. We have assumed a distance to the LMC and, more specifically, to SN 1987A of 51.4 ± 1.2 kpc (Panagia et al. 1991, later updated in Panagia 1999), whereas the photometric properties of the instrument were taken from the WFPC2 Instrument Handbook (Heyer & Biretta 2004; namely, inverse sensitivity PHOTFLAM =1.461 × 10⁻¹⁶ and photometric zero-point ZP(Vegamag) =17.564). We derive a median value of the luminosity of the 189 objects with Hα excess of 4 × 10³¹ erg s⁻¹ or 10⁻² L_☉. A histogram of the individual L(Hα) values is shown in Figure 5 (dashed line).

Zoom In Zoom Out Reset image size

It should be noted that, because of the width of the Hα filter, ΔHα includes a small contribution due to [N ii] emission features at 6548 Å and 6584 Å. Generally, their intensities do not exceed 1.2% and 3.5%, respectively, of the Hα line intensity, and in many instances the [N ii] doublet lines are much weaker, i.e., <0.1% of the Hα line intensity (e.g., Edwards et al. 1987; Hartigan et al. 1995). These intensities should be even lower for PMS stars in the LMC because of the lower metallicity relative to the Milky Way. However, to take a conservative approach, we have simply adopted an average value that is half of the maximum for Galactic T-Tauri stars and we give as uncertainty half of the measured range. We will therefore adopt average values of (1.8 ± 1.8)% and (0.6 ± 0.6)% relative to the Hα intensity, for the 6584 Å and 6548 Å [N ii] lines, respectively (note that the uncertainties after the "±" sign here are meant to show the entire span of the range, not just the 1σ value). Considering the width of the Hα filter of the WFPC2, only the 6548 Å [N ii] line contaminates the Hα measurements, whereas for measurements made with the Hα filter on-board the ACS, both lines will be included (albeit only in part for the 6548 Å line). Although the resulting corrections on the Hα intensity are quite small (namely, 0.994 for the WFPC2 F656N filter and 0.979 for the ACS F658N filter, on average), they are systematic and it is therefore good practice to apply them in all cases.

Typically, the combined total uncertainty on our L(Hα) measurements is ∼15% and is completely dominated by the statistical uncertainty on the Hα photometry. The systematic uncertainty on the distance to the field of SN 1987A accounts for 5% (see Panagia 1999), whereas the typical uncertainty on the absolute sensitivity of the instrumental setup is of order 3%.

3.2. The Hα Equivalent Width

If the stars defining the reference template had no Hα absorption features (we have established above that they do not have Hα in emission), their V − Hα index would correspond to that of the pure continuum and could then be used to derive the equivalent width of the Hα emission line, W_eq(Hα). While this approximation might be valid for stars with a conspicuous Hα emission, it is clearly not applicable in general. It is, however, possible in all cases to derive the level of the Hα continuum by using properly validated models atmospheres.

As discussed above, the models of Bessell et al. (1998) reproduce quite reliably the observed broadband colors (see, e.g., Figure 2). This means that, even though the models might lack the resolution necessary to realistically reproduce spectral lines (hence the small discrepancy between the dashed line and the squares in Figure 3), the level of the continuum can be trusted. One can, therefore, fit the continuum in those models in a range containing the Hα line (e.g., 6500–6620 Å) and fold the resulting "Hα-less" model spectra through the instrumental setup (e.g., with Synphot), so as to derive the magnitude of the sole continuum in the Hα band (Hα^c) as a function of spectral type or effective temperature. We present in the Appendix the relationships between V, I, and Hα^c for several HST instruments and filters. The difference between the observed Hα magnitude and Hα^c provides a direct measure of W_eq(Hα), as we explain below.

We recall that the equivalent width of a line is defined as

$$\begin{equation} W_{\rm eq}= \int {(1 - P_\lambda) \, d\lambda }, \end{equation} \tag{ 3 }$$

where P_λ is the profile of the line, or the spectrum of the source divided by the intensity of the continuum, and the integration is extended over the spectral region corresponding to the line. If the line profile is very narrow when compared to the width of the filter, i.e., if the line falls completely within the filter bandpass as is typical of emission lines, then W_eq(Hα) is simply given by the relationship:

$$\begin{equation} W_{\rm eq} ({\rm H}\alpha) = \mathrm{RW} \times [1-10^{-0.4 \times ({\rm H\alpha -H\alpha }^c)}], \end{equation} \tag{ 4 }$$

where RW is the rectangular width of the filter (similar in definition to the equivalent width of a line), which depends on the characteristics of the filter. Values of RW for the Hα bands of current and future HST instruments are given in the Appendix (Table 4). As in the case of L(Hα), the statistical uncertainty on the equivalent width is dominated by the uncertainty on the Hα magnitude, typically of order 15% in our selection. Moreover, the value of W_eq(Hα) obtained in this way is also subject to some systematic uncertainties produced by the model atmospheres used to interpolate the level of the continuum, such as, e.g., metallicity, surface gravity, or spectral resolution, and should therefore be taken with caution. On the other hand, we show in the Appendix that the equivalent widths (in absorption) obtained by applying this method to the high-resolution model spectra of Munari et al. (2005) are in excellent agreement with those obtained via standard spectrophotometry (see Figure 18).

Similarly, the values of W_eq(Hα) that we obtain in the field of SN 1987A are fully consistent with those of young PMS stars. In Figure 6, we show the value of W_eq(Hα) as a function of V − I for the 4156 selected stars. As mentioned earlier, there are 189 objects with a ΔHα excess above the 4σ level and they are indicated as thick symbols (circles and triangles) in Figure 6. Since the equivalent width shown in Figure 6 is that of the pure emission component, the spectra of stars with small W_eq, say ≲10 Å, actually show an Hα in absorption. For this reason, we conservatively ignore stars with measured W_eq(Hα) < 10 Å, since this value is about the largest absorption equivalent width expected for normal stars (see Figure 18 in the Appendix). A total of 164 objects satisfy this condition, but, in light of our conservative threshold, this must be considered as a lower limit to the number of objects with genuine Hα excess.

Zoom In Zoom Out Reset image size

Values of W_eq(Hα) for the sample range from 10 Å to 650 Å, with a median of 19 Å. These values are in excellent agreement with those obtained by Panagia et al. (2000) from the R − Hα color excess of the stars in this field and represent a refinement of those measurements, since these are based on a more accurate determination of the level of the continuum. These equivalent widths are typical of PMS stars (e.g., Muzerolle et al. 1998b) and about 2 orders of magnitude smaller than those expected for pure nebular emission (see discussion in Section 6). While objects bluer than V − I ≃ 0.2 are probably Ae/Be stars, there are 133 redder stars with significant Hα excess (marked as circles in Figure 6) that are most likely bona fide low-mass PMS stars of various ages, as we will show in the following sections. A histogram of their Hα luminosities is shown in Figure 5 (solid line).

It is generally accepted (see, e.g., Hartmann et al. 1998) that the energy L_acc released by the magnetospheric accretion process goes toward ionizing and heating the circumstellar gas. In this hypothesis, the bolometric accretion luminosity L_acc can be determined from a measurement of the reradiated energy and, in particular, via the measured Hα luminosity that is produced in the process.

In order to determine the exact relationship between L(Hα) and L_acc, we have used literature values of L_acc and L(Hα) measurements as recently summarized by Dahm (2008) for PMS stars in the Taurus–Auriga association. As pointed out by Dahm himself, there may be a serious problem in determining an empirical relationship between L_acc and L(Hα) measurements because the two quantities were determined from non-simultaneous observations. The intensity of the Hα line from PMS sequence stars is known to vary by about 20% on a timescale of a few hours and by as much as a factor of 2–3 in a few days (e.g., Fernandez et al. 1995; Smith et al. 1999; Alencar et al. 2001). For this reason, a relationship between L(Hα) and L_acc is necessarily rather uncertain, although a very clear trend can be seen in the L_acc versus L(Hα) plot shown by Dahm (2008). His logarithmic best fit would provide a slope of 1.18 ± 0.26 for such a relationship. On the other hand, theoretical models (e.g., Muzerolle et al. 1998b) would predict logarithmic slopes of about unity for low accretion rates, hence for faint Hα luminosities, and shallower slopes for higher luminosities.

In the absence of compelling evidence in favor or against a steep slope, we will adopt a logarithmic slope of unity for the empirical L_acc versus L(Hα) relationship, i.e., a constant ratio L_acc/L(Hα), and determine the proportionality constant from an elementary fit of the data summarized by Dahm (2008). On this basis we obtain

$$\begin{equation} \log L_{\rm acc} = (1.72 \pm 0.47) + \log L({\rm H}\alpha). \end{equation} \tag{ 5 }$$

As anticipated, the resulting uncertainty on L_acc is rather large, as high as a factor of 3, but this is the best one can obtain with the available data sets. In view of the convenience and the feasibility of determining L_acc from L(Hα) for PMS outside our own Galaxy, it would be of paramount importance to obtain a proper calibration of the L_acc versus L(Hα) relationship by making simultaneous spectral observations of PMS stars over a suitably wide wavelength baseline, i.e., covering at least the range 3000–8000 Å, and possibly extending it to the near-IR with measurements of Paschen and Brackett lines. High quality data of this type are already available for the Orion Nebula Cluster, collected by dedicated projects both from ground-based observatories and with the HST (i.e., the Legacy Program on Orion; P.I.: Massimo Robberto) and their analysis is in progress to provide a more accurate calibration of the L_acc versus L(Hα) relationship (N. Panagia et al. 2010, in preparation).

The median value of the accretion luminosity thus obtained for our sample of 133 PMS stars is ∼0.57 L_☉, in good agreement with the value of 0.54 L_☉ derived from the U-band excess measured by Romaniello et al. (2004) for stars in the same field, especially considering that their selection of stars with excess emission is less stringent than our 4σ level (we discuss this point in Section 7).

Once L_acc is known, the mass accretion rate follows from the free-fall equation that links the luminosity released in the impact of the accretion flow with the rate of mass accretion $\dot{M}_{\rm acc}$ according to the relationship:

$$\begin{equation} L_{\rm acc} \simeq \frac{{\it GM}_*\,\dot{M}_{\rm acc}}{R_*} \left(1 - \frac{R_*}{R_{\rm in}}\right), \end{equation} \tag{ 6 }$$

where G is the gravitational constant, M_* and R_* are the mass and photospheric radius of the star, and R_in is the inner radius of the accretion disk. The value of R_in is rather uncertain and depends on how exactly the accretion disk is coupled with the magnetic field of the star. Following Gullbring et al. (1998), we adopt R_in = 5 R_* for all PMS objects.

As regards R_*, we derive it from the luminosity and effective temperature of the stars, which in turn come from the observed color (V − I) and magnitude (V), properly corrected for interstellar extinction as provided by Romaniello (1998), Panagia et al. (2000), and Romaniello et al. (2002). The adopted distance to SN 1987A is 51.4 kpc (Panagia et al. 1991; Panagia 1999), as mentioned above. The stellar mass M_* was estimated by comparing the location of a star in the H–R diagram with PMS evolutionary tracks. As for the latter, we adopted those of Degl'Innocenti et al. (2008), for metallicity Z = 0.007 or about one-third Z_☉, as appropriate for the LMC (e.g., Hill et al. 1995; Geha et al. 1998). These new PMS tracks were specifically computed with an updated version of the FRANEC evolutionary code (see Chieffi & Straniero 1989 and Degl'Innocenti et al. 2008 for details and Cignoni et al. 2009 for an application). However, in order to assess how differences in the evolutionary models affect our results, in Section 6.1 we also consider PMS tracks from other authors, namely, those of D'Antona & Mazzitelli (1997) and Siess et al. (2000).

Combining Equations (5) and (6), we can link the mass accretion rate $\dot{M}_{\rm acc}$, in units of M_☉ yr⁻¹, to L(Hα):

$$\begin{eqnarray} \log \frac{\dot{M}_{\rm acc}}{M_\odot \, {\rm yr}^{-1}} & = &-7.39 + \log \frac{L_{\rm acc}}{L_\odot } + \log \frac{R_*}{R_\odot } - \log \frac{M_*}{M_\odot } \nonumber\\ & = & (-5.67 \pm 0.47) + \log \frac{L({\rm H}\alpha)}{L_\odot } \nonumber\\ && + \log \frac{R_*}{R_\odot } - \log \frac{M_*}{M_\odot }. \end{eqnarray} \tag{ 7 }$$

The distribution of mass accretion rates found in this way for the 133 PMS stars in our sample is shown in Figure 7. The median value of 2.9 × 10⁻⁸ M_☉ yr⁻¹ is in remarkably good agreement with the median rate of 2.6 × 10⁻⁸ M_☉ yr⁻¹ found by Romaniello et al. (2004) from the U-band excess.

Zoom In Zoom Out Reset image size

The uncertainty on $\dot{M}_{\rm acc}$ is dominated by the approximate knowledge of the ratio of L_acc and L(Hα) in Equation (5). As explained above, empirical measurements and theoretical models suggest a value of L_acc/L(Hα) ≃ 52 but with an uncertainty of a factor of 3. Since the ratio, however, is the same for all stars, the comparison between different objects is not hampered by this uncertainty, as long as the statistical errors are small.

As for the other quantities in Equation (7), we discuss here the sources of statistical uncertainty, while systematic errors are addressed separately in Section 6. With our selection criteria, the typical uncertainty on L(Hα) is 15% and is dominated by random errors. The uncertainty on R_* is typically 7%, including a 5% systematic uncertainty on the distance modulus. As for the mass M_*, since it is determined by comparing the location of a star in the H–R diagram with evolutionary tracks, both systematic and statistical uncertainties are important. The uncertainty on the temperature is mostly statistical and of order 3%, while that on the luminosity is 7%, comprising both random errors (1% uncertainty on the bolometric correction and 3% on the photometry) and systematic effects (5% uncertainty on the distance modulus). When we interpolate through the PMS evolutionary tracks to estimate the mass, the uncertainties on T_eff and L imply an error of ∼7% on M_*. However, an even larger source of systematic uncertainty on M_* comes from the evolutionary tracks, as we explain in Section 6. In summary, the combined statistical uncertainty on $\dot{M}_{\rm acc}$ is 17%.

Having identified a population of bona fide PMS stars, with well defined accretion luminosities and mass accretion rates, it is instructive to place them in the H–R diagram. We do so in Figure 8 where bona fide PMS stars are marked as large dots, with a size proportional to their $\dot{M}_{\rm acc}$ value, according to the legend. Also shown in the figure are the theoretical isochrones from the FRANEC models of Degl'Innocenti et al. (2008) for Z = 0.007 and ages as indicated (in units of Myr, from the stellar birth line; see Palla & Stahler 1993). The dashed line, corresponding to a 40 Myr isochrone, defines in practice the zero-age MS.

Zoom In Zoom Out Reset image size

It is noteworthy that the majority of PMS objects in Figure 8 are rather close to the MS and would have easily been missed if no information on their Hα excess had been available. Indeed, it is customary to identify PMS stars in a color–magnitude diagram by searching for objects located above and to the right of the MS, since this is where one would expect to find very young objects. This method was first used by Gilmozzi et al. (1994) to identify a population of PMS stars outside the Milky Way, namely, in the LMC cluster NGC 1850. More recent applications of this method include, e.g., those of Sirianni et al. (2000), Nota et al. (2006), and Gouliermis et al. (2007).

Unfortunately, this method of identification is not very reliable, since the existence of an age spread and the presence in the same field of a considerably older population, as well as unaccounted patchy absorption, all tend to fill up the parameter space between the MS and the birth line, thereby thwarting any attempts to identify PMS stars on the basis of their effective temperature and luminosity alone. In the specific case of the SN 1987A field studied here, Panagia et al. (2000) have shown that the age spread is remarkable, with several generations of young stars with ages between 1 and 150 Myr superposed on a much older field population (0.6–6 Gyr). In such circumstances, as Figure 8 shows, broadband photometry alone would have not identified as such most of our bona fide PMS stars. In accordance with their proximity to the MS, the majority of bona fide PMS stars turn out to be relatively old. The ages of individual objects were determined by interpolating between the isochrones in the H–R diagram. They range from 1 to ∼50 Myr, with a median age of 13.5 Myr, in very good agreement with the estimated age of both SN 1987A and its nearby Star 2 of ⩽13 Myr (see Scuderi et al. 1996, and references therein). Panagia et al. (2000) find that about 35% of the stars younger than 100 Myr in this field have a typical age of 12 Myr, again consistent with the estimated age of the bona fide PMS stars.

The relatively old age of our PMS stars might appear at odds with the rather large mass accretion rates that we obtain for that population, namely, ∼10⁻⁸ M_☉ yr⁻¹. According to current models of viscous disk evolution (see Hartmann et al. 1998; Muzerolle et al. 2000), the mass accretion rate of a ∼14 Myr old PMS population should be of order 3 × 10⁻¹⁰ M_☉ yr⁻¹. Those models, shown in Figure 9 as a solid line, adequately reproduce the trend of decreasing $\dot{M}_{\rm acc}$ with stellar age observed for very low mass Galactic T-Tauri stars (see, e.g., Calvet et al. 2004), but underpredict our measurements (indicated by diamonds) by over an order of magnitude.

Zoom In Zoom Out Reset image size

On the other hand, recent observations of PMS stars in the cluster Trumpler 37 by Sicilia-Aguilar et al. (2006), from whom the data points in Figure 9 are taken, indicate that the value of $\dot{M}_{\rm acc}$ for ∼10–20 Myr old G-type PMS stars (large circles) is considerably larger, indeed of order ∼10⁻⁸ M_☉ yr, in excellent agreement with our results. With a median mass of 1.2 M_☉ and a median V − I color of 0.56, our bona fide PMS stars are fully consistent with the G spectral type and suggest that stars of higher mass have a higher $\dot{M}_{\rm acc}$ value at all ages. As mentioned in Section 1, Muzerolle et al. (2003, 2005), Natta et al. (2004, 2006), and Calvet et al. (2004) suggest a dependence of the type $\dot{M}_{\rm acc} \propto M^2$ from the analysis of a sample of low-mass and intermediate-mass T-Tauri stars, albeit with a large uncertainty on the index (see also Clarke & Pringle 2006). A similar trend is proposed by Sicilia-Aguilar et al. (2006), yet with even lower statistical significance due to the limited mass range covered by their observations. The range of masses that we cover here is also quite limited and does not allow us to address this issue in detail. We will, however, do so in a forthcoming paper (G. De Marchi et al. 2010, in preparation), where we compare the accretion rates in this field with those of a much younger and an order of magnitude more numerous population of PMS stars in the Small Magellanic Cloud. While in the SN 1987A field the detection limits imposed by the WFPC2 photometry do not allow us to reach stars with mass below ∼0.9 M_☉, in that paper, based on ACS observations, we will study stars down to ∼0.4 M_☉.

The main sources of systematic uncertainty on the derived mass accretion rates are (1) discrepancies in the isochrones and evolutionary tracks, (2) reddening, (3) Hα emission generated by sources other than the accretion process, and (4) the contribution of the nebular continuum to the colors of the stars. We discuss all of these effects in this section.

6.1. Models of Stellar Evolution

6.2. Reddening

Dust extinction systematically displaces stars in the H–R diagram toward lower luminosities and effective temperatures. This is shown graphically in Figure10, for six model stars, three with T_eff = 6000 K (large circles) and three with T_eff = 7000 K (large squares). The effects of a reddening increase from E(V − I) = 0 to E(V − I) = 1 on the location of the stars in the H–R diagram is shown by the smaller symbols, corresponding to increments of E(V − I) = 0.1. We have adopted here the extinction law as measured specifically in the field of SN 1987A by Scuderi et al. (1996).

Zoom In Zoom Out Reset image size

Since the reddening vector in the H–R plane is almost parallel to the lines of constant radius (as shown by the dot-dashed line for R = 10 R_☉), extinction affects only marginally the estimate of the stellar radius. However, at low temperatures the reddening vector crosses the evolutionary tracks (solid lines), thus leading to an underestimate of the stellar mass, if no extinction correction is applied. In turn, this will result in a systematically higher value of the R_*/M_* ratio and, therefore, in an overestimate of the mass accretion rate derived via Equation (6).

In order to better gauge the uncertainty on the value of R_*/M_* caused by unaccounted extinction, we have determined the mass and radius of each one of the artificially reddened model stars shown in Figure 10, in exactly the same way in which we did that for bona fide PMS stars, and compared those with the known input values. The result is shown graphically in Figure 11, where, the relation between the derived and input R_*/M_* ratios is shown as a function of color excess E(V − I).

Zoom In Zoom Out Reset image size

Regardless of the input parameters (L, T_eff, or R_*/M_*), it appears that over the range explored here the derived R_*/M_* ratio increases slowly with color excess (see dashed line that best fits the median values indicated by the crosses). In particular, underestimating the E(V − I) color excess by ∼0.4 mag would lead to a 30% overestimate of R_*/M_* (and hence of $\dot{M}_{\rm acc}$). This appears unlikely, considering that the best estimate of the total extinction—including the Milky Way component—toward Star 2, one of the two bright close companions to SN 1987A, is E(V − I) = 0.21 (Scuderi et al. 1996), with a 1σ dispersion over the field of ∼0.07 (Romaniello 1998; Panagia et al. 2000). According to Zaritsky et al. (2004), the average reddening toward cool (5500 K <T_eff < 6500 K) stars in the LMC is A_V = 0.43, thus implying E(V − I) ≃ E(B − V) ≃ 0.14. Therefore, for an average LMC star-forming region, omitting the extinction correction would result in a 10% overestimate of the mass accretion rate. This error, albeit systematic, is smaller than most other systematic uncertainties and comparable to the measurement errors.

Nonetheless, we would like to stress here that the effects of patchy absorption can be more severe in regions of high extinction. Therefore, whenever possible, one should apply extinction corrections to each individual stars, as it was done for the photometric catalog of Romaniello et al. (2002) that we use here.

6.3. Hα Emission not Due to the Accretion Process

The underlying assumption in the determination of the mass accretion rate of PMS stars is that the gravitational energy L_acc released in the accretion process goes into heating of the gas at the boundary layer (see, e.g., Hartmann et al. 1998). In this case, the luminosity of the Hα emission line can provide a measure of the accretion energy, because it acts as a natural "detector" of the luminosity released in the accretion process.

Thus, the issue to address is whether the measured L(Hα) is produced only and entirely as a result of the energy released in the accretion process. Obviously, depending on the geometry and orientation of the circumstellar disk and on the specifics of the accretion process, some of the ionizing energy could escape without effectively ionizing the local gas. This would result in an underestimate of the true mass accretion rate for some stars, but this is true of any method that uses H emission lines as indicators of the accretion process and, therefore, it does not affect the comparison of our results with those of others (see Figure 9).

One possibility that we can readily exclude is that the Hα emission that we detect is due to chromospheric activity. The total chromospheric emission of a solar-type star on the subgiant branch corresponds to a luminosity of about 10⁻⁴ L_☉ (Ulmschneider 1979; Pasquini et al. 2000), i.e., 2 orders of magnitude less than what we measure.

On the other hand, we must also consider the possibility that Hα emission might occur in discrete H knots along the line of sight, in places unrelated with the PMS object. Furthermore, hot nearby massive stars may ionize the H gas in which PMS stars are immersed, eventually producing Hα emission when the gas recombines. Both these effects would result in a overestimate of $\dot{M}_{\rm acc}$, but we show below that the probability of this happening is very small.

As for intervening Hα emission along the line of sight, it would have to arise in knots of ionized H, e.g., very compact H ii regions, contained within our photometric apertures. However, if such knots exist, they should not project only against stars and, therefore, narrowband images should reveal Hα emission where broadband imaging shows no or very low signal. Since this is not observed, statistical arguments rule out this possibility as a significant source of the observed Hα excess emission.

While it is possible that some of the Hα emission that we detect is due to diffuse nebular emission in the H ii region not powered by the accretion process, if the emission is extended and uniform over an area of at least 05 radius around the star its contribution cancels out with the rest of the background when we perform aperture photometry (using an aperture of 02 radius for the object and an annulus of 03–05 for the background). The subtraction would not work if the emission were not uniform, as, for example, in the case of a filament that projects over the star but that does not cover completely the background annulus. However, we have visually inspected all objects with excess Hα emission and excluded the few cases where a spurious filament of this type could be seen (see Romaniello 1998; Panagia et al. 2000).

Another case in which the emission does not cancel out is when the gas ionized by an external source is associated to the star and completely contained within our photometric aperture. We recall here that the typical size of a circumstellar disk is of order 100 AU (e.g., Hartmann et al. 1998), or about 0.02 pixel at the distance of SN 1987A, and that the median L(Hα) that we measure for our PMS stars is 4 × 10³¹ erg s⁻¹, or about 10⁻² L_☉ (see Figure 5). We can then calculate how close a PMS stars should be to a hot, UV-bright star for the latter to flood the PMS circumstellar disk with enough ionizing radiation to cause at least 10% of the Hα luminosity that we measure from that source.

Let r = 100 AU be the external radius of the circumstellar disk and d the distance to an ionizing star that illuminates the disk, assumed to have zero inclination (i.e., face-on), so that ϕ = r²/(4 d²) is the solid angle subtended by the disk. Let β₂ be the recombination coefficient to the first excited level in an optically thick H gas. If V is the volume of the gas and n_e and n_p are the electron and proton densities per unit volume, the global ionization equilibrium equation (case B; see Baker & Menzel 1938) is then:

$$\begin{equation} N_{L} \, \phi = \beta _2 \, n_e \, n_p \, V, \end{equation} \tag{ 8 }$$

where N_L is the rate of ionizing photons produced by the external star. Denoting with α_3,2 the Hα emission coefficient for the same gas, the luminosity of the Hα recombination line can be written as

$$\begin{equation} L({\rm H}\alpha) = \alpha _{3,2} \, n_e \, n_p \, V, \end{equation} \tag{ 9 }$$

which allows us to express the luminosity of the Hα emission caused by the external star as a function of its ionizing photon rate N_L:

$$\begin{equation} L^{\rm ext}({\rm H}\alpha) = \alpha _{3,2} / \beta _2 \, N_L \,\phi. \end{equation} \tag{ 10 }$$

The minimum distance beyond which the contribution of the ionizing star to the total Hα emission becomes negligible (i.e., 10% or less) is given by

$$\begin{equation} d^2_{\rm min} = \frac{r^2}{4} \, \frac{\alpha _{3,2}}{\beta _2} \, \frac{10 \, N_L}{L({\rm H}\alpha)}. \end{equation} \tag{ 11 }$$

Table 1 gives the value of ϕ for a typical Hα luminosity of 10⁻² L_☉, as a function of the spectral type of the ionizing star. Here, we have used the compilation of Panagia (1973) for the appropriate stellar parameters. The table also provides the minimum distance d_min in pc and in arcsec (at the distance of SN 1987A) for a circumstellar disk of radius r = 100 AU. We recall here that 1 arcsec corresponds to 10 WFPC2 pixels. The brightest star in the field studied here is an object with log T_eff = 4.7 and log(L/L_☉) = 5.8, with an implied mass of ∼60 M_☉ and an age less than 1 Myr (see Panagia et al. 2000). The corresponding value of N_L, from Panagia (1973), is 4 × 10⁴⁹ s⁻¹, thus implying d_min = 0.95 pc or 38 WFPC2 pixel for a gas temperature of 10,000 K. A careful inspection of the images shows that none of our candidate bona fide PMS stars falls closer to this star than d_min. The same applies to the second and third brightest stars in the field, respectively, with log T_eff = 4.7, log(L/L_☉) = 5 and log T_eff = 4.67, log(L/L_☉) = 4.75, as well as to all fainter objects.

Table 1. Possible Effects of Nearby Ionizing Stars

SpT	Teff	log(L/L☉)	log NL	ϕ	dmin
	(K)		(s−1)		(pc)	('')
O4	50,000	6.11	49.93	3.42e-08	1.31	5.26
O6	42,000	5.40	49.08	2.43e-07	0.49	1.97
O7	36,500	4.81	48.62	6.93e-07	0.29	1.17
O9	34,500	4.66	48.08	2.43e-06	0.16	0.62
B0	30,900	4.40	47.36	1.27e-05	0.07	0.27
B1	22,600	3.72	45.29	1.50e-03	0.01	0.03

Notes. The values of ϕ and d_min are given for a 100 AU circumstellar H disk with a temperature of 10,000 K. For gas temperatures of 5000 K and 20,000 K, the value of ϕ is, respectively, higher and lower by 7.5%, while d_min is 3.6% larger and smaller, respectively.

Download table as:  ASCIITypeset image

It is convenient to calculate, for each PMS object, the total Hα luminosity L^ext(Hα) caused by all neighboring ionizing stars, expressed as

$$\begin{equation} L^{\rm ext}({\rm H}\alpha) = \frac{r^2}{4} \, \frac{\alpha _{3,2}}{\beta _2} \, \sum _{i=1}^N \frac{(N_L)_i}{d_i^2}, \end{equation} \tag{ 12 }$$

where (N_L)_i is the ionizing photon rate of star i and d_i is its distance from the PMS object being considered. Besides the three brightest stars mentioned above, we have included as potential sources of the detected Hα emission an additional 12 stars in the field with log T_eff ⩾ 4.4 and log(L/L_☉) ⩾ 4, as classified by Romaniello et al. (2002). We note that these estimates of L^ext(Hα) are firm upper limits because we are assuming that (1) all disks are face-on relative to the hot stars, and that (2) the projected distances on the sky are indeed the actual distances from the hot stars. The values of L^ext(Hα) obtained in this way for the 133 bona fide PMS stars span the range 10⁻⁶–10⁻⁴ L_☉, or typically 3 orders of magnitude less than the median PMS Hα luminosity, as shown by the histogram in Figure 12. We can, therefore, conclude that in this specific field the contribution of hot ionizing stars to the measured excess Hα emission is negligible. We must stress, however, that the effect might be much more pronounced in younger star-forming regions, such as the Orion Nebula or 30 Dor.

Zoom In Zoom Out Reset image size

6.4. Nebular Continuum

We have developed and successfully tested a new self-consistent method to reliably identify all PMS objects actively undergoing mass accretion in a resolved stellar population, regardless of their age, without requiring spectroscopy. The method combines broadband V and I photometry with narrowband Hα imaging to (1) identify all stars with excess Hα emission using as a reference template of the photospheric luminosity in the Hα band the average V − Hα color of stars with small photometric errors, (2) convert the excess Hα magnitude into Hα luminosity L(Hα) via the absolute photometric calibration of the Hα band, (3) convert, if desired, the excess Hα magnitude into the Hα emission equivalent width, using model spectra to fit the level of the continuum in the Hα band, (4) derive the accretion luminosity L_acc from L(Hα) using a relationship based on recent literature values, and (5) obtain the mass accretion rate $\dot{M}_{\rm acc}$ from the free-fall equation that links it to L_acc via the stellar parameters (mass and radius) that we derive from the H–R diagram using PMS tracks for the appropriate metallicity.

Since no spectroscopy of individual objects is needed, observations requiring just a few hours of exposure time can suffice to reveal hundreds of PMS stars simultaneously, determine their Hα luminosity with an uncertainty of less than 15%, and derive their mass accretion rates with the same accuracy as allowed by conventional spectral line analysis.

As regards the Hα luminosity, the virtue of this new method is that it derives the reference luminosity of the stellar photosphere inside the specific Hα band simply from the average V − Hα color of stars with the same V − I index (see Figure 4). This is possible because, at any given time, the majority of the stars in a stellar field have no excess Hα emission. Obviously, the average luminosity of the stellar photosphere in the Hα band reflects more or less pronounced Hα absorption features, depending on the spectral type of the objects, and precisely for this reason it defines the reference luminosity with respect to which any excess emission should be sought and measured. The uncertainty associated with the derived L(Hα) is dominated by statistical errors in the photometry and does not exceed 15% for our 4σ selection (see Section 3.1). Note that this approach is also applicable to populations comprising mostly very young objects, such as the Orion Nebula, since excess Hα emission is found in only about one-third of the PMS stars at any given time and, therefore, the majority of the very young objects will not have a conspicuous Hα emission, thus correctly defining the reference template (see Section 3.1).

In spite of the good accuracy on L(Hα), deriving the mass accretion rate requires the use of a rather uncertain relationship linking L(Hα) to the accretion luminosity (see Section 4 and also Dahm 2008). The uncertainty on the relationship can have a systematic effect of up to a factor of 3 on the value of $\dot{M}_{\rm acc}$. This is the unfortunate consequence of our incomplete understanding of the physics of the accretion process, plagued by the lack of direct measurements of the bolometric accretion luminosity L_acc, as well as by the paucity of simultaneous measurements of the hot continuum excess emission and of the emission line fluxes in Hα or other lines (e.g., Paβ and Brγ). Although annoying, the systematic uncertainty of up to a factor of 3 on $\dot{M}_{\rm acc}$ is similar to that plaguing age determinations based on the comparison with theoretical isochrones, which can reach a factor of 2 (see Section 6.1).

On the other hand, these systematics affect in the same way all determinations of $\dot{M}_{\rm acc}$ that make use of emission lines diagnostics, including those based on a detailed study of the profile and intensity of the Hα line (e.g., Muzerolle et al. 1998a), since their calibration relies on the same measurements. For this reason, the relative values of $\dot{M}_{\rm acc}$ that we obtain are trustworthy and the comparison of our results with those obtained with other methods remains meaningful.

As discussed in Sections 3 and 4, as an application of our method, we have studied the field of SN 1987A, where we identify 133 bona fide PMS stars, based on their large Hα excess. The median value of the mass accretion rate that we find for these stars, namely, 2.9 × 10⁻⁸ M_☉ yr⁻¹, is in very good agreement with the median value of 2.4 × 10⁻⁸ M_☉ yr⁻¹ found by Romaniello et al. (2004) from the analysis of the distribution of the (U − B) color excess in the same field. While we have identified 133 bona fide PMS stars with an Hα excess above the 4σ level, Romaniello et al. conclude that 765 stars with an excess U-band luminosity larger than ∼0.035 L_☉ could be bona fide PMS candidates, albeit at the 1.5σ level.

The approach followed by Romaniello et al. (2004) is statistical in nature and, while providing a solid estimate of the average mass accretion rate, it does not allow one to reliably identify individual PMS stars. Therefore, Romaniello et al. turned to the Hα excess emission (namely, R − Hα) as a more direct indicator of ongoing accretion and found a strong correlation between the U-band continuum excess and W_eq(Hα) for stars with U − B>0.5 and W_eq(Hα)>3 Å. Not surprisingly, the more stringent threshold on the equivalent width adopted in the present study, i.e., W_eq(Hα)>10 Å, results in the identification of a smaller number of bona fide PMS stars, but also in a much smaller uncertainty. Indeed, the method that we present here is more direct than the one based on the U-band excess, since it does not require the use of model atmospheres to predict the level of the intrinsic photospheric continuum in the U band from the observed optical colors, and as such it results in a more reliable determination of the accretion luminosity.

Our method, therefore, makes it possible for the first time to carry out a quantitative and systematic study of the properties of a large number of PMS stars, including those in dense and/or distant star-forming regions. This allows us to tackle some still open questions as to the evolution of both the stars and their disks. For instance, in the field of SN 1987A, we have studied how the mass accretion rate changes with time as our bona fide PMS stars approach their MS. We find that $\dot{M}_{\rm acc}$ decreases more slowly with time than what is observed for low-mass T-Tauri stars in the Galaxy (see Section 5) and predicted by models of viscous disk evolution (Hartmann et al. 1998). However, our results are in excellent agreement with the (so far limited) measurements of the mass accretion rate of G-type stars in the Milky Way (e.g., Sicilia-Aguilar et al. 2006), thus lending credit to the hypothesis that $\dot{M}_{\rm acc}$ might depend on the stellar mass.

Moreover, our method allows us to investigate whether and how the accretion process is affected by the chemical composition and density of the molecular clouds or by the proximity of hot, massive ionizing stars. We will address the effects of metallicity in a forthcoming paper (G. De Marchi et al. 2010, in preparation) devoted to a comparison of star formation in the Small Magellanic Cloud and LMC. As for the effect of hot massive stars, we have studied the distribution of the 133 bona fide PMS objects with respect to the baricenter of the 15 hottest stars in the field, described in Section 6.3. We find a clear anti-correlation, shown in Figure 13, between the frequency of the identified PMS stars and their distance from the ionizing objects (marked with blue star symbols in Figure 14). This effect is not seen for non-PMS objects of similar brightness, whose distribution is uniform over the whole field, confirming that the trend is not due to problems in detecting faint objects near the brightest stars. The trend remains the same whether the baricenter of the hot stars is computed using as weight their bolometric luminosity or their ionizing flux. This anti-correlation is quite remarkable in several ways, as we discuss below.

Zoom In Zoom Out Reset image size

In Figure 14, we show a schematic map of the distribution of all 133 PMS objects with respect to the 15 ionizing stars. The baricenter of the ionizing flux is shown by the large square and concentric annuli are drawn from it, in increments of 35''. For reference, a red ellipse marks the position of SN 1987A. Symbols of different colors indicate PMS stars with different Hα luminosity, from red for L(Hα)>2 × 10⁻² L_☉ to white for <8 × 10⁻³ L_☉, with green diamonds corresponding to intermediate values. An inspection to Figure 14 immediately reveals two important facts: (1) there are very few PMS stars near the baricenter of the massive objects and (2) their Hα luminosity is systematically lower than that of PMS stars farther away.

Panagia et al. (2000) had already pointed out that there is no overdensity of young PMS objects around hot massive stars in this field. This shows that low-mass and massive stars do not follow each other, suggesting that their formation mechanisms are separate and distinct. But what Figure 14 additionally shows is that there are very few PMS stars with strong Hα luminosity (red squares) around the young massive stars, with only 20% of them within a 70 arcsec radius of the baricenter. Conversely, stars with low Hα luminosity (white squares) are distributed rather uniformly across the field.

Zoom In Zoom Out Reset image size

With a median age of 13.5 Myr, most PMS objects are much older than the 1–2 Myr old hot stars and must belong to a previous generation. Thus, no spatial relationship (whether a correlation or anti-correlation) should be expected between the distribution of the two types of objects. On the other hand, while we cannot exclude that originally fewer PMS stars formed in the region now populated by young massive stars, there are indications that the local paucity of PMS objects is only apparent and caused by the presence of the hot young stars.

The fact that PMS objects near the youngest hot stars are both less numerous and fainter in Hα emission suggests that their circumstellar disks have been considerably eroded and made less efficient by enhanced photo-evaporation. We have shown earlier (see Section 6.3 and Figure 12) that the instantaneous contribution to the Hα luminosity from the ionizing radiation of nearby massive stars is negligible. The results shown in Figure 14 confirms that this is the case. If it were not and the observed Hα luminosity were mostly due to fluorescence of the circumstellar disks under the effect of the ionizing radiation of massive stars, the radial distribution should be opposite to the one observed, with brighter Hα luminosities for PMS stars closer to the baricenter of the youngest hot stars. Thus, if the apparent paucity of PMS objects in their vicinity is the result of photo-evaporation, the latter must be mostly due to the effects of non-ionizing far-ultraviolet (FUV) radiation. When integrated over the lifetime of massive stars (∼2 Myr), their FUV radiation can be quite significant for the disruption of nearby circumstellar disks. In particular, this also means that there might be many more objects still in their PMS phase, close to the hot stars, but we have no way to identify them as such, since the photo-evaporation of their disks has stopped the accretion process.

Recent theoretical studies of the erosion of circumstellar disks via photo-evaporation due to external FUV radiation (e.g., Adams et al. 2004) or to the PMS star itself (e.g., Alexander et al. 2006a, 2006b) suggest timescales of order 8–10 Myr for the complete dissolution of the disks. Although this is the right order of magnitude for the timescale of the process that we see, there are still some discrepancies. For instance, it is clear that if disks were to completely dissolve under the effect of their own central star within 5–10 Myr, as suggested by Alexander et al. (2006b), the majority of our PMS objects, with a median age of 13.5 Myr, should not show any sign of accretion. Similarly, with a typical age of 2 Myr, the hottest stars in our field would have not had time to deposit enough UV photons on to the PMS disks to seriously erode them, if the process is to take of order 10 Myr (Adams et al. 2004). While recently Clarke (2007) suggested that the effect of nearby massive young stars may be stronger, leading to disk exhaustion timescales of order 2 Myr, her calculations apply to the immediate vicinity (<1 pc) of very massive stars, such as the O6 star θ₁C in Orion. The density of massive stars in the field of SN 1987A is much lower than that present in the Orion Nebula and the distances to our PMS objects typically an order of magnitude larger than those considered by Clarke (2007).

Nevertheless, while it is clear that the theory of disk photo-evaporation does not yet capture all the details, we believe that it can greatly benefit from the observations presented here. In particular, it will be necessary to explain which mechanisms allow isolated circumstellar disks to feed their central PMS stars for well over 20 Myr, at least in a low-metallicity, low-density environment such as the field of the LMC. At the same time, the models should justify the quick (≲3 Myr) demise of the PMS disks within a (projected) radius of ∼10 pc of a handful of late O stars. Extending our investigation, as we plan to do, to more and denser regions, such as NGC 1850 and 30 Dor in the LMC and NGC 346 and NGC 602 in the Small Magellanic Cloud, will offer increased statistics and a much wider variety of star formation environments.

We are very thankful to Pier Giorgio Prada Moroni and Scilla Degl'Innocenti for sharing with us their pre-main-sequence evolutionary tracks ahead of publication and to Massimo Robberto for useful discussions. G.D.M. is grateful to STScI and ESO for the hospitality, via their Science Visitor programs, during the preparation of this paper. We are also indebted to an anonymous referee for precious comments and suggestions that have helped us to improve the presentation of our work.

Over the years, a large number of observations of stellar populations have been collected with the HST/WFPC2 and HST/ACS through the V-like filters F555W and F606W and the I-like filter F814W. We have therefore computed color relationships between the magnitude in the R-like band F675W (m₆₇₅) and the V − I-like colors m₅₅₅ − m₈₁₄ and m₆₀₆ − m₈₁₄. To this aim, we have used both observed and theoretical stellar spectra, taken, respectively, from the BPGS Spectrophotometry Atlas (see Gunn & Stryker 1983) and from the stellar atmosphere models of Bessell et al. (1998).

The BPGS atlas contains nearby bright dwarfs and giants of near-solar metallicity, whereas from the library of Bessell et al. (1998) we have selected specifically those models with metallicity index [M/H] = −0.5, as appropriate for the LMC (Dufour 1984), surface gravity g = 4.5 and effective temperatures in the range 3500  K ⩽T_eff⩽ 40,000 K. Using the standard HST synthetic photometry package Synphot (Laidler et al. 2008), we have folded these spectra through the instrumental response of the WFPC2, the ACS/WFC, the ACS/HRC, and the WFC3/UVIS cameras, in order to obtain the expected magnitudes in the bands listed above. For comparison with ground-based photometry, we have also calculated the expected magnitudes through the Johnsn–Cousin V, R, and I bands.

We show in Figure 15, for the WFPC2 filter system, the predicted m₅₅₅ − m₆₇₅ color as a function of m₅₅₅ − m₈₁₄ for all the stars in the BPGS atlas. The distribution is remarkably linear and the dashed line represents the best fit in the color range 0 < m₅₅₅ − m₈₁₄ < 1.5 (filled symbols), corresponding to the relationship m₆₇₅ = 0.34 × m₅₅₅ + 0.66 × m₈₁₄. The color–color diagram in the same bands based on the models of Bessell et al. (1998) is shown in Figure 16 and the relationship is practically indistinguishable, namely, m₆₇₅ = 0.35 × m₅₅₅ + 0.65 × m₈₁₄. In this case, the dashed line shows the least-squares fit to stars in the range 4000 K ⩽T_eff⩽ 10,000 K, marked by the filled symbols. As already noted in Section 2, it appears that neither the luminosity class nor the metallicity has a noticeable impact on the color relationship in these bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This is also the case when we consider the m₆₀₆ − m₈₁₄ color, for which the relationship with m₆₇₅ is given in Table 3. In the same table, we also list the results for the ACS/WFC, ACS/HRC, and WFC3/UVIS cameras, as well as those for the standard Johnson–Cousin system normally used from the ground (note that the R and I bands in this system do not coincide with the R and I bands of Johnson's original system, as mentioned in Section 2). As expected, the coefficients in the color relationship vary considerably with the band, and differences exist also between the same band on various HST cameras (or, at least, bands with the same name).

As briefly discussed in Section 3, the magnitude corresponding to the level of the continuum in the Hα band can be derived by means of properly validated model atmospheres, such as those of Bessell et al. (1998) that we have shown to reproduce quite accurately the observed broadband colors (see Figure 2). This means that, even though the models might lack the spectral resolution necessary to realistically reproduce spectral lines (hence the small discrepancy between the dashed lines and the squares in Figure 2), the level of the continuum can be trusted. We have, therefore, fitted the continuum in the Bessell et al. (1998) spectra over the range 6500–6620 Å and have replaced the actual spectrum in that range with the interpolated continuum. Finally, we have used Synphot to fold the model spectra, modified in this way, through the instrumental setup and derive in this way the continuum Hα magnitudes (Hα^c) for the various Hα filters on-board the HST, as a function of the spectral type or T_eff of the star.

With this information it is now possible to build color–color diagrams like those of Figures 15 and 16, and look for a color relationship between Hα^c and V − I. In Figure 17, we show the specific case of the WFPC2. The m₅₅₅ − m^c₆₅₆ color varies very smoothly and very slowly with the effective temperature and is a remarkably linear function of the m₅₅₅ − m₈₁₄ color over the whole range, and particularly for stars with 4000 K <T_eff < 10,000 K (filled circles) with which we are mostly concerned in this work. Similar relationships exist for other bands and photometric systems and they are given in Table 3. We note that the constant 0.351 in the color relationship accounts for the conspicuous Hα absorption line in the spectrum of Vega, which defines the standard for the VEGAMAG magnitude system.

Zoom In Zoom Out Reset image size

As indicated in Section 3, the difference between the observed Hα magnitude and Hα^c can be converted to the equivalent width E_eq(Hα) by means of Equation (4) and using the values of RW as listed in Table 4. To prove this, we have applied Equation (4) to the model atmospheres of Munari et al. (2005), which provide spectra with a resolving power of 11,500 and a 1 Å sampling, and compared the results to literature measurements of Hα equivalent widths (in absorption).

Ducati (1981) compiled an extensive homogeneous catalog of Hα line measurements, from photographic and photoelectric techniques, for about 2300 stars. We have converted Ducati's index α to W_eq(Hα) following the calibration provided by Strauss & Ducati (1981) and Cester et al. (1977). The values of W_eq(Hα) obtained in this way are shown in Figure 18, where crosses correspond to dwarfs and boxes to giants. Ducati's (1981) stars are classified according to their spectral type, which we have converted to V − I following Bessell (1990).

Zoom In Zoom Out Reset image size

As for the application of Equation (4) to the models of Munari et al. (2005), this was done for model spectra with 4000 < T_eff/K < 40,000, log g = 4.5, and solar metallicity and for the WFPC2 photometric system, namely, F555W for V, F814W for I and F656N for Hα. It should be noted that, for stars with broad Hα lines (around spectral type A), the equivalent width determined in this way is systematically smaller than the true one, since part of the line wings fall outside of the Hα filter. However, the measurements of Ducati (1981) are subject to exactly the same effect, since they were derived in a similar fashion and with a Hα filter of comparable effective width. The value of W_eq(Hα) obtained in this way is shown in Figure 18 (solid line), as a function of V − I. The agreement with the observations is very satisfactory.