Protoplanetary Disks around Sun-like Stars Appear to Live Longer When the Metallicity is Low^*

In their final formation stages, after dissipating most of the envelope of dust and gas in which they were born, Sun-like stars still have a mostly gaseous disk from which they continue to accrete while the star contracts until hydrogen ignites in the core (C. Hayashi 1966; F. H. Shu et al. 1987; F. Palla & S. W. Stahler 1999; R. B. Larson 2003). It is in these disks that planetary systems form (P. J. Armitage 2010). In nearby Galactic star-forming regions with solar composition, the disks are short-lived, with decay timescales of ∼2 Myr and an accretion phase that rapidly tapers off (K. E. Haisch et al. 2001; J. Hernandez et al. 2008; E. E. Mamajek 2009; D. Fedele et al. 2010; A. J. W. Richert et al. 2018).

Since the majority of stars in the Universe formed at cosmic noon around redshift z ∼ 2 (P. Madau et al. 1996; S. Lilly et al. 1996), it is important to investigate the lifetime of protoplanetary disks in the lower-metallicity conditions (Z ∼ 0.1 Z_⊙) characteristic of those environments (e.g., M. Zwaan et al. 2005; J. Fynbo et al. 2006). In this respect, it has been suggested that lifetimes in lower-metallicity environments should be even shorter than in the solar neighborhood because disks are less shielded by dust from the X-ray emission of the central star (B. Ercolano et al. 2008, 2009). The fraction of disk-bearing stars with near-infrared (NIR) excess measured in two small low-metallicity (∼0.1 Z_⊙) star-forming regions in the extreme outer Galaxy initially (C. Yasui et al. 2009, 2010) appeared to support this hypothesis (but see Section 6).

In contrast, photometric studies of equally low-metallicity star-forming regions in the Magellanic Clouds using the Hubble Space Telescope (HST) have revealed thousands of stars that appear to systematically accrete at higher rates and for times exceeding 30 Myr (G. De Marchi et al. 2010, 2011a, 2013, 2017; K. Biazzo et al. 2019; S. Tsilia et al. 2023). In particular, observations (E. Sabbi et al. 2007) of the massive star-forming region NGC 346, in the Small Magellanic Cloud (SMC), revealed approximately ∼690 objects with strong Hα excess emission and Hα equivalent widths (EWs) larger than 20 Å (G. De Marchi et al. 2011a). As such, these objects are candidate pre-main-sequence (PMS) stars actively accreting matter from their circumstellar disks. G. De Marchi et al. (2011a) find that these PMS objects form two distinct populations, with a clearly bimodal age distribution: about half of the objects are younger than ∼2 Myr and the other half older than ∼8 Myr, including over 200 objects with ages in the range of 15−30 Myr (see Section 5 for details).

Since in Galactic star-forming regions accreting PMS stars older than 15 Myr are very rare (G. J. Herczeg & L. A. Hillenbrand 2014), it has been proposed (G. De Marchi et al. 2010, 2011a, 2013, 2017) that the PMS accretion phase lasts longer in the low-metallicity environments typical for the Magellanic Clouds (∼0.1−0.3 Z_⊙), and by inference at cosmic noon. While other processes, such as chromospheric activity or binary interactions, might be responsible for the Hα excess emission in HST photometry of stars hotter than ∼7500 K, this is unlikely because the measured Hα emission EW is more than an order of magnitude higher than for chromospherically active binaries (D. Montes et al. 1995). However, without a means to obtain spectra for these objects, the nature of the Hα excess emission of the sources in NGC 346 has so far remained unclear.

With the NIRSpec Micro Shutter Array (MSA), we studied a sample of these candidate PMS objects in NGC 346, as well as normal main-sequence (MS) stars with no detectable Hα excess emission. Here, we present the spectra and discuss the physical properties of 10 representative sources with similar masses (∼0.9−1.8 M_⊙), T_eff in the range of 4500−8000 K, and spanning a range of ages (∼0.1−30 Myr).

The structure of this paper is as follows. The observations are presented in Section 2, while the extraction and calibration of the spectra are discussed in Section 3. The spectral analysis is contained in Section 4, where we derive the physical parameters of the objects. A detailed discussion of the age of PMS stars and their uncertainties is provided in Section 5. Conclusions are presented in Section 6.

The central regions of NGC 346 were observed with NIRSpec on 2022 July 26 using the G235M/F170LP medium-resolution grating (R ∼ 1000) covering the range of 1.7−3.2 μm (proposal number 1227). A total of 40 sources were observed simultaneously through the NIRSpec MSA, of which half had been identified as candidate PMS stars (G. De Marchi et al. 2011a), owing to their Hα excess emission, and the other half as candidate MS stars (no detectable excess in Hα). The targets were selected from a catalog of about 80,000 sources in this field detected and measured with HST photometry (E. Sabbi et al. 2007), out of which about 690 have Hα excess emission exceeding 5σ (G. De Marchi et al. 2011a).

Using the eMPT software tool (N. Bonaventura et al. 2023), we identified the most suitable targets (and corresponding MSA configuration) for the specific orientation of the telescope at the time of the observations. The MSA contains about 250,000 microshutters, each one defining a slit with a width of 020 and a length of 046, while the pitch between the microshutters is 027 in the spectral direction and 053 in the spatial direction (P. Ferruit et al. 2022). We elected to reserve for each source a "slitlet" composed of three neighboring microshutters in the cross-dispersion direction, in order to also properly characterize the nebular background around our targets. Each slitlet spans a region of 02 × 15 on the sky. Several constraints were applied to derive the final list of 40 targets, one of which was to avoid sources with contaminants within three microshutters for any position of the source in the slitlet. Indicating with m_I the I-band magnitude of the possible target source in the catalog by E. Sabbi et al. (2007), we considered as viable targets only those for which any contaminants within the three microshutters were fainter than m_I + 5. We further constrained the eMPT tool to exclude objects for which the wavelengths corresponding to the Paα or Brβ spectral lines would fall in the gap between detectors to ensure that both lines can be measured simultaneously. For some sources, the wavelength range around the Brγ line is also recorded.

Each source was observed at three different nodding positions in the slitlet by slewing the telescope in the cross-dispersion direction by an amount corresponding to the microshutter pitch (053). Each target star was observed for 481 s in each of the three microshutters, at the same relative position, for a total 1443 s of integration. Due to the fixed MSA pattern, it is not possible to locate all targets at the nominal center of a microshutter.

In this first work, we present the results for 10 representative sources with similar masses (0.9–1.8 M_⊙, median 1.2 M_⊙) that span a range of ages, including both younger and older PMS candidate stars, as well as MS objects. The stars were selected to have a signal-to-noise ratio (SNR) >3 (as measured per pixel with the DER_SNR routine; F. Stoehr et al. 2008), to be free from contaminants in the immediate vicinity of the selected microshutters (some of the contaminants were not detectable in the optical HST images but are visible in the NIRCam exposures; O. C. Jones et al. 2023; N. Habel et al. 2024), and for which also the Brγ is recorded in the spectrum (only for one of the objects this is not the case). In Table 1, we list the coordinates and key stellar parameters. The masses and ages of the PMS candidates were determined by G. De Marchi et al. (2011a) by comparing the effective temperature T_eff and bolometric luminosity L_bol of the stars, derived from extinction-corrected V and I broadband HST photometry, to theoretical PMS evolutionary tracks (S. Degl'Innocenti et al. 2008; E. Tognelli et al. 2011) for metallicity Z = 1/8 Z_⊙, as appropriate for NGC 346 (for more details, see G. De Marchi et al. 2011a and Section 5). For candidate MS stars, which were selected among objects with no Hα excess emission in the catalog of G. De Marchi et al. (2011a), the masses were derived in the same way, by comparing the T_eff and L_bol of the stars to MS isochrones from the same set of models, after reddening correction. (Note that older PMS candidates and MS stars generally share the same position in the color–magnitude diagram (CMD), the only difference being the presence or absence of Hα excess emission with EW >20 Å.) The resulting masses are in the range of ∼1.2−1.3 M_⊙, comparable with those of the PMS candidates.

Table 1. Properties of the Stars in Our Sample

ID	R.A.	Decl.	Age	M*	AV	$\mathrm{log}\,T$	$\mathrm{log}\,L$	R*	m200	WPaα	WBrβ	WBrγ	$\mathrm{log}\,{L}_{{\rm{a}}}^{{\rm{H}}\alpha }$	$\mathrm{log}\,{L}_{{\rm{a}}}^{\mathrm{Br}\gamma }$	$\mathrm{log}\,\dot{M}$
			(Myr)	(M⊙)	(mag)	(K)	(L⊙)	(R⊙)	(mag)	(Å)	(Å)	(Å)	(L⊙)	(L⊙)	(M⊙ yr−1)
9884	14.8247272	−72.1544850	0.1	0.9	0.55	3.65	1.10	5.8	20.4	−142.1	−41.8	−21.5	1.02	1.28	−5.30
9685	14.8206427	−72.1862390	1.5	1.8	0.60	3.74	0.96	3.3	21.6	−65.1	−27.1	−10.6	0.25	0.39	−6.74
24897	14.8067222	−72.1913565	5.5	1.2	0.65	3.75	0.42	1.7	23.2	−210.0	−302.0	−25.8	−0.09	−0.27	−7.51
12559	14.8044053	−72.1929698	12.5	1.3	0.66	3.91	0.76	1.3	23.3	−157.8	−389.3	−8.0	−0.27	−0.99	−8.38
20408	14.7178504	−72.1690186	16.0	1.2	0.56	3.89	0.47	1.0	24.0	−165.2	−139.6	...	−0.32	...	...
26651	14.8120579	−72.1876153	20.0	1.0	0.60	3.82	0.31	1.1	23.9	−195.0	−129.7	−7.6	−0.59	−1.41	−8.76
40571	14.8088971	−72.1823255	27.0	0.9	0.55	3.80	0.03	0.8	24.4	−387.4	−290.0	−45.4	−0.10	−0.41	−7.68
10147	14.8854390	−72.1667293	...	1.3	0.56	3.89	0.86	1.5	23.1	...	...	...	...	...	...
14090	14.8811044	−72.1616790	...	1.4	0.55	3.88	0.68	1.3	23.4	...	...	...	...	...	...
19639	14.8660216	−72.1617760	...	1.2	0.42	3.84	0.45	1.2	23.9	...	...	...	...	...	...

Note. The columns correspond to the source identification (ID), J2000.0 coordinates (R.A., decl.), age, stellar mass (M_*), extinction (A_V ), effective temperature ($\mathrm{log}\,T$), bolometric luminosity ($\mathrm{log}\,L$), stellar radius (R_*), magnitude in the F200W NIRCam band (m₂₀₀), EW of the Paα, Brβ, and Brγ lines, accretion luminosity derived from the L(Hα) and L(Brγ) luminosity, and mass accretion rate ($\mathrm{log}\,\dot{M}$). Stellar parameters are from G. De Marchi et al. (2011a), NIRCam F200W AB magnitudes are from N. Habel et al. (2024; F277W AB magnitudes for stars 9685 and 40571). All spectral lines have negative EW, indicating emission.

Download table as:  ASCIITypeset image

The exposures were processed with the NIRSpec ramp-to-slope pipeline, developed by the ESA NIRSpec Science Operations Team (SOT), which performs the following basic data reduction steps: bias subtraction, reference pixel subtraction, linearity correction, dark subtraction, and finally count-rate estimation, including jump detection and cosmic-ray rejection (S. M. Birkmann et al. 2022). From the count-rate images, the flux- and wavelength-calibrated spectra were obtained using Stage 2 of the SOT NIRSpec Instrument Pipeline Software (NIPS; C. Alves de Oliveira et al. 2018b) to perform the following default operations (C. Alves de Oliveira et al. 2018a): subtract background, extract subimage containing the spectral trace, assign wavelength and spatial coordinates to each pixel therein, flat-fielding, and finally flux calibration. This procedure generates a rectified spectrum resampled on a regular two-dimensional (2D) grid, which in the case of our filter and grating combination (F170LP/G235M), produces images that span 7 × 1414 pixel², covering the wavelength range of 1.66−3.16 μm, and 046 in the cross-dispersion direction. The 2D-rectified spectra of two of our sources are shown in Figure 1.

Zoom In Zoom Out Reset image size

To optimize the SNR, we computed the one-dimensional (1D) spectrum of each star from the 2D-rectified product by coadding three pixel rows around the nominal source position, as given by the MSA planning tool. To account for any slit losses and to perform absolute calibration, we rescaled each spectrum in order to match the corresponding magnitude of each source in the F200W band as observed with NIRCam (N. Habel et al. 2024). The scale factor was derived by folding each observed spectrum through the NIRCam F200W band using the stsynphot tool (STScI Development Team 2020) and comparing the resulting synthetic magnitude with the one actually observed. For stars 9685 and 40571, the magnitude in the F277W band was used instead since these objects were not observed with the F200W filter. The NIRCam magnitudes are listed in Table 1.

Given the highly variable nebulosity in the region, background subtraction is a critical step. We take as a reference the two strong He i lines at 1.87 μm and 2.06 μm, which are fully nebular in nature (the mass and T_eff of our targets are such that these lines are not expected to originate in their photospheres). We measure the flux of the He i lines by integrating over their width, after subtraction of the continuum. This is done for the spectrum containing source+background and for the two background spectra sampled by the microshutters just above and below the source. We take as a template of the nebular background the spectrum from the microshutter where the integrated He i line flux is closest to that in the spectrum with the source. We then rescale the background template (by multiplication) until the integrated He i line flux agrees with that in the spectrum of the source, and finally subtract the rescaled template from the spectrum of the source. The multiplicative factor is always smaller than 1.2 (i.e., 20%) and typically ranges between 1.0 and 1.1, indicating small intensity variations in the nebular background around our sources, which our procedure corrects.

In Figure 1, the examples of sources 9884 and 20408 illustrate the two cases of uniform and highly varying backgrounds. The outline of the three microshutters in the slitlet is shown on top of images in different bands, namely, Hα (F658N) from the HST Advanced Camera for Surveys observations (E. Sabbi et al. 2007) and F444W from JWST NIRCam observations (O. C. Jones et al. 2023; N. Habel et al. 2024). Since there are three nodding positions, in total five microshutters are shown in projection.

Our final goal is to correct the H recombination lines in our sources, i.e., to remove the signatures of the nebular background. Therefore, the uncertainty introduced by this procedure depends on how well the observed He i intensity variations match the variations in the H recombination lines. To assess this, we compare to one another the spectra of the two background microshutters (those above and below the source, namely, "b" and "d" in Figure 1), after having registered them to have the same He i line flux. For the three prominent H recombination lines in our study, intensity variations range from 1% to 10%, with a typical difference of 5% for Paα, 6% for Brγ, and 5% for Brβ. For some targets, the He i line flux in the background spectra matches that in the source spectrum (within Poisson uncertainties), so no rescaling of the background spectrum is needed before subtraction. This is the case when the nebular background is rather constant around the source, like for source 9884, where the background level in the Hα band is rather uniform in the four regions probed by the nodding pattern around the source. The situation is very different for source 20408, where simple background subtraction without rescaling would not be appropriate.

We finally note that, by design, our nodding strategy places the spectrum of the source and the background spectra at the same location on the detector and in the field of view, thereby avoiding uncertainties caused by possible flat-fielding and throughput mismatches as well as nonresponsive pixels.

After background subtraction, the spectra extracted from the three nodding positions of each source were median averaged to improve the SNR, and pixels flagged as saturated, nonlinear, or affected by high detector noise were excluded.

The background-subtracted and combined spectra are shown in Figure 2, while the derived physical parameters of the stars are listed in Table 1. The typical statistical uncertainty on the extracted spectra, after all pipeline steps and background subtraction, is ∼3 × 10⁻²⁰ erg cm⁻² s⁻¹ Å⁻¹. This is shown by the shaded areas around each spectrum in Figure 2. The corresponding SNR of the spectra (as defined by F. Stoehr et al. 2008) ranges from 3 for source 40571 to 42 for source 9884. The typical uncertainty (median value) and range (from the 17th to 83rd percentiles) of the measured line luminosities and EW are 2% (2%−5%) for Paα, 6% (4%−12%) for Brβ, and 9% (2%−16%) for Brγ and are dominated by the uncertainty on the level of the underlying continuum. The absolute calibration of the spectra is based on NIRCam photometry in the F200W or F277W bands, which for our sources in the range of 20 < m₂₀₀ < 24 has a typical uncertainty of 0.1 mag or 10% (N. Habel et al. 2024).

Zoom In Zoom Out Reset image size

All PMS candidates show prominent H recombination lines in their spectra, in particular Paα and Brβ (Brγ is also in emission). For some of the PMS candidates, Brδ and Br are also detected. Conversely, the spectra of MS stars show H absorption lines but no emission lines. No He i nebular lines are present in any of these spectra, confirming good nebular background subtraction (see Section 3). In Figure 3, we compare the spectrum of a PMS candidate of intermediate age (source 12559, ∼12 Myr, black line) with a representative MS spectrum (red line) obtained as the mean of the three MS stars in our sample. Three H absorption lines are clearly present in the MS spectrum (Br, Paα, Brδ), in contrast with the prominent emission features of the PMS star.

Zoom In Zoom Out Reset image size

In this section, we discuss the physical properties derived from the spectra, starting from those of MS stars. We also discuss and quantify three types of spectral features that are characteristic of young stellar objects, namely, atomic H recombination lines, molecular H excitation lines, and continuum NIR excess.

4.1. Metallicity, Effective Temperature, and Surface Gravity

4.2. H Recombination Lines and Mass Accretion

4.3. H₂ Excitation Lines and Outflows

The background-subtracted spectra of the PMS stars in Figure 2 reveal a conspicuous number of molecular hydrogen (H₂) excitation lines. In particular, the lines corresponding to the rovibrational transitions v = 1 → 0 S(0), 1 → 0 S(1), 1 → 0 S(2), 1 → 0 S(3), 2 → 1 S(1), 2 → 1 S(2), 2 → 1 S(3), as well as 1 → 0 Q(1) through to 1 → 0 Q(7) are observed. Since the spectra are background subtracted, the molecular gas must be spatially associated with the stars.

By knowing the exact location of the stars in the microshutter, we can set constraints on the distribution of the H₂ molecular material. Figure 4 shows the spectra of source 40571, without background subtraction; the side panel indicates the position of the source inside the microshutter (as provided by the MSA configuration file), which is intermediate between pixel rows 2 and 3 in the extracted rectified trace (note that rows 1 and 7 of the trace are only partly illuminated). While the He i and Brγ nebular emission lines extend across the full height of the slit, the spatial extent of the v = 1 → 0 S(1) emission line is more limited around the specific location of the source. The same is true for the other six PMS candidates, indicating that H₂ is likely associated with circumstellar material.

Zoom In Zoom Out Reset image size

The extent of the H₂ emission in the spatial direction is about 01 (1 pixel) and is not resolved. Therefore, we can only set an upper limit to the extent of the H₂ circumstellar structure (possibly an outflow in the disk or a filament), which at the distance of NGC 346 corresponds to ∼6000 au. This is clearly just an upper limit, and we cannot set any constraints on the shape or structure at the origin of the H₂ emission, but we note that this size is within a few times the typical ≲1000 au size of disks around Galactic PMS stars measured in Taurus (J. R. Najita & E. A. Bergin 2018).

Some of these H₂ lines are observable from the ground in nearby protostars and their intensity ratios are used to characterize the physical conditions and excitation mechanism in their circumstellar environments (T. P. Greene et al. 2010). Here, we provide a preliminary characterization of the environments and excitation mechanisms from the intensity ratios of the S(1) lines in the v = 2 → 1 and v = 1 → 0 transitions. For star 12559, the intensity ratio is 0.47, suggesting nonthermal ultraviolet electron pumping, for which R. Gredel & A. Dalgarno (1995) predict a ratio of 0.54 (the value is higher if extinction is present). The intensity ratio of the S(3) lines in the v = 2 → 1 and v = 1 → 0 transitions agree with this scenario. For stars 9685, 24897, and 40571, the ratios are in the range of 0.11−0.18, characteristic of LTE molecular gas at temperatures of ∼2200−2700 K, as commonly seen in shocks from outflows. Stars 26651, 9984, and 20408 were analyzed for the 2 → 1 S(3) line ratio because the S(1) line was affected by bad pixels or not covered in the wavelength range of the spectra. These three stars also have line ratios consistent with the shock excitation of H₂ gas in LTE at ∼2200−2700 K.

In summary, the intensity ratios of the S(1) and S(3) lines in the v = 2 → 1 and v = 1 → 0 transitions suggest that the likely excitation mechanism for one of the stars (12559) is ultraviolet electron pumping, while for the other objects shocks appear to be the cause. The measured H₂ line ratios in these stars point to excitation mechanisms that likely trace a variation in underlying physical processes from outflows and circumstellar disk material commonly seen in young stars in the Milky Way.

4.4. NIR Excess and Circumstellar Dust

Finally, the spectral energy distributions (SEDs) of the PMS stars in our sample present NIR excess caused by dust emission, which provides additional evidence of the presence of circumstellar disks. In Figure 5, we compare the SEDs of the stars obtained by fitting the continuum of the spectra in Figure 2 using a linear interpolation with a step of ∼0.1 μm. The derived SEDs were corrected for reddening, using the A_V values in Table 1 with the extinction law of G. H. Rieke & M. J. Lebofsky (1985), and are normalized at 1.75 μm for comparison. The gray dashed line corresponds to the average MS star spectrum shown in Figure 5 and closely matches the model spectra (M. S. Bessell et al. 1998) for metallicity [M/H] = −0.5, surface gravity $\mathrm{log}g=4.5$ and T_eff of 4500 K, 6500 K, and 8000 K (see the legend). This comparison confirms that in this wavelength and temperature range the spectral shape in the Planck tail of the distribution is not particularly sensitive to the exact T_eff value. The differences observed in the SEDs of the PMS stars, with respect to those of the MS objects, clearly indicate the presence of NIR excess. This further confirms that the stars with excess emission in the H recombination lines are not standard MS stars. If their continuum matched that of MS stars, the emission lines might simply be the result of an incomplete subtraction of the nebular background. Instead, the markedly different continuum shape provides independent evidence that these objects are not MS stars and are still surrounded by circumstellar material.

Zoom In Zoom Out Reset image size

All PMS candidates except for source 12559 reveal NIR excess also in the NIRCam photometric bands F115W, F200W, and F277W (N. Habel et al. 2024), with excess values ranging from 0.25 to 1.3 mag compared to the colors of the MS stars in our sample. In general terms (see Table 1), younger PMS candidates have typically larger NIR excess, with some scatter. The scatter, however, is not surprising: although the mass of the circumstellar disk decreases over time, the amount of NIR excess does not directly scale with age because of the effects introduced by different disk geometry and morphology, including inclination (A. Garufi et al. 2018). For instance, the presence of spirals and shadows is associated with a high NIR excess, while the presence of rings in the disks appears to cause low NIR excess (A. Garufi et al. 2018), irrespective of the actual PMS age. Therefore, a pronounced NIR excess for an older PMS candidate, like star 20408, can be expected.

In this section, we discuss the ages of the PMS stars in our sample, how they were derived from the photometry and their uncertainties. Some elements of this discussion have been presented by G. De Marchi et al. (2011a), but here we address the matter in more detail and through a number of independent arguments, all of which point to an old age for a significant fraction of the accreting PMS stars in NGC 346.

A search for PMS stars in NGC 346 (G. De Marchi et al. 2011a) revealed approximately ∼690 candidates with strong Hα excess emission in their photometry, corresponding to Hα EW >20 Å. In the V, V − I CMD, after reddening correction, these PMS candidates form two distinct groups, a red group and a blue group, with a median V − I color separation of 0.75 mag at all magnitudes. When the T_eff and L_bol of these stars, derived from the photometry, are compared in the Hertzsprung–Russell (H-R) diagram with widely used PMS models (L. Siess et al. 2000; E. Tognelli et al. 2011; J. Tang et al. 2014), the 0.75 mag color separation translates into a clearly bimodal age distribution, with about half of the objects younger than ∼2 Myr and the other half older than ∼8 Myr, including over 200 stars with ages in the range of 15−30 Myr. To help readers appreciate the systematically different locations of the younger and older PMS candidates in the CMD and H-R diagram, we have included in the left and right panels of Figure 6, respectively, the CMD and the H-R diagram as published by G. De Marchi et al. (2011a), with those authors' permission. These figures are relevant because they were used to select the spectroscopic targets for the investigation reported here.

Zoom In Zoom Out Reset image size

Hereafter, we will refer to these two groups of PMS stars as "older PMS candidates" and "younger PMS candidates," respectively. The former include the stars with Hα excess emission in Figure 6 with color 0.1 ≲ V − I ≲ 0.8, which G. De Marchi et al. (2011a) showed to be older than ∼8 Myr (median age 20 Myr), to have masses in the range of 0.7−1.5 M_⊙, and to be located primarily in the outskirts of the association. The latter include objects with Hα excess emission in Figure 6 with color 0.8 ≲ V − I ≲ 2.0, younger than ∼2 Myr (median age 0.9 Myr), with masses in the range of 0.4−4.0 M_⊙, and located near the center of the association, in compact subclusters.

In the following subsections, we address three specific points, which are however intertwined, namely, the location of the PMS sources in the CMD, their reddening, and their spatial distribution. We also explore how representative the older emitting PMS sources are of the older PMS population in NGC 346 as a whole.

5.1. CMD and Photometric Uncertainties

Here, we discuss the location of the PMS sources in the CMD and show that their small photometric uncertainties, combined with the low metallicity of the NGC 346 cluster, allow us to set stringent constraints on the relative ages of the younger and older PMS stars.

The typical photometric uncertainty in the V and I bands for the PMS stars in the catalog of E. Sabbi et al. (2007) used by G. De Marchi et al. (2011a) is 0.02 mag, which translates into typical uncertainties of about 1%−2% on T_eff and of ∼6% on L_bol. In that work, the value of T_eff was derived by comparing the observed V − I color, after correction for reddening, with those stemming from the model atmosphere of M. S. Bessell et al. (1998) for metallicity [M/H] = −1.0 and gravity $\mathrm{log}g=4.5$, as appropriate for the PMS sources in NGC 346, and for the specific HST bands used in those observations (see G. De Marchi et al. 2010, 2011a for details). The stellar parameters derived spectroscopically in Section 4 attest to the accuracy of the effective temperature obtained from the photometry. Bolometric luminosities are based on the adopted SMC distance modulus of 18.92 ± 0.03, corresponding to ∼61 kpc (R. W. Hilditch et al. 2005; S. C. Keller & P. R. Wood 2006). While it has long been known that the SMC extends about 5 kpc along the line of sight (e.g., L. T. Gardiner & M. R. S. Hawkins 1991; D. Nidever et al. 2011; P. Yanchulova Merica-Jones et al. 2017; C. E. Murray et al. 2024), there is no evidence that this is the case for the specific NGC 346 star-forming region. The diameter of the region is about 60 pc, and one expects the extent along the line of sight to be similar. This will have no measurable effects on the dispersion in the photometry and hence on the derived value of L_bol.

Figure 6 shows that older PMS candidates with Hα excess emission (which here we show to have excess emission also in the Paschen and Brackett lines) appear very close to the location of the MS, shown by the solid gray line in the right panel of Figure 6. Objects of this type are not observed in nearby Galactic star-forming regions. However, their position in the H-R diagram is fully consistent with theoretical models of stellar evolution in low-metallicity environments. At the metallicity of NGC 346 (Z = 1/8 Z_⊙ or Z ≃ 0.002), theoretical evolutionary tracks (S. Degl'Innocenti et al. 2008; E. Tognelli et al. 2011) show that a 1 M_⊙ star reaches the MS after about 20 Myr, as we illustrate in Figure 7. The time is about 27 Myr for a 0.9 M_⊙ object like source 40571 in Table 1. At solar metallicity, models by the same authors show that this time is twice as long, or about 43 Myr for a 1 M_⊙ star. We verified that both the MIST (A. Dotter 2016) and PARSEC (A. Bressan et al. 2012) evolutionary tracks would produce very similar results.

Zoom In Zoom Out Reset image size

More importantly, even if the lifetime of circumstellar disks did not depend on metallicity, Figure 7 implies that in low-metallicity environments there is a higher probability of finding PMS stars close to the MS, which are still actively accreting. Indeed, photometric investigations of star-forming regions across the Magellanic Clouds have systematically shown this result over the past decade (G. De Marchi et al. 2010, 2011a, 2011b, 2013, 2017; K. Biazzo et al. 2019; R. Carini et al. 2022; S. Tsilia et al. 2023; M. Vlasblom & G. De Marchi 2023). Obviously, if disk lifetimes increase at lower metallicity, as our observations appear to indicate (see Section 6), the likelihood of finding still accreting PMS stars close to the MS is even higher.

Figure 7 presents an additional important implication for uncertainties on ages, which are derived through isochrone comparison and as such are affected by uncertainties on T_eff and L_bol. The latter in turn depend on uncertainties in the photometry and on the extinction correction (typical uncertainties on $\mathrm{log}{T}_{\mathrm{eff}}$ and $\mathrm{log}{L}_{\mathrm{bol}}$ are shown by the error bars in Figure 7). None of these uncertainties depend on metallicity. Since at lower metallicity stars transit more quickly across the H-R diagram, a given uncertainty on $\mathrm{log}{T}_{\mathrm{eff}}$ and $\mathrm{log}{L}_{\mathrm{bol}}$ will translate into a lower uncertainty on the age. We quantify this in Figure 7, where PMS evolutionary tracks for 1 M_⊙ stars and two metallicities are compared. At Z = 1/8 Z_⊙, the star leaves the bottom of the Hayashi track at t = 4 Myr and reaches the MS region 12 Myr later, having become about 1500 K warmer, at an average rate of 150 K Myr⁻¹. At solar metallicity, the time spent on the radiative part of the PMS track is 24 Myr and T_eff increases by about 1300 K, so the average transfer rate is about 55 K Myr⁻¹. During the transfer phase, which accounts for the majority of the PMS lifetime, age uncertainties are about three times smaller at Z = 1/8 Z_⊙ than at solar metallicity (the actual factors, for a 1 M_⊙ star, are 2.7, 2.5, and 2.6 for the PISA, MIST, and PARSEC tracks, respectively). This analysis shows that PMS stars in NGC 346 form two distinct age groups.

5.2. Reddening

5.3. Spatial Distribution of Younger and Older PMS Stars

5.4. Accreting and Nonaccreting Older PMS Stars

Thanks to the excellent combination of sensitivity and spatial resolution of NIRSpec on JWST, we have been able to observe spectroscopically a sample of extragalactic solar-mass candidate PMS stars for the first time. The candidates were identified through Hα excess emission in their photometry. Our new data set confirms that these objects are indeed PMS stars, with the spectral signatures of active accretion and the presence of colder molecular gas and of dust in circumstellar material. The spectroscopic confirmation of the PMS nature of the photometrically selected candidates suggests that it will be possible to extend the use of narrowband imaging methods to reliably identify and study star-forming objects in a variety of environments, also in galaxies further away in the Local Group.

The spectral features presented in Section 4 are common to both younger and older PMS stars in our sample, with ages spanning from ∼0.1 to ∼30 Myr. Independently of the age spread, all our objects reveal typical PMS signatures: prominent H recombination lines implying large accretion luminosities, H₂ excitation lines suggesting molecular material in the immediate vicinity of the stars, and NIR excess witnessing the presence of dusty circumstellar disks. This suggests that in NGC 346 circumstellar disks live longer than typical disks in the Milky Way and PMS stars can accrete at sustained rates for over 20 Myr. This is in contrast with what is observed in the Milky Way, where at solar metallicity, disks around stars of ∼1 M_⊙ dissipate very quickly (S. Pfalzner et al. 2022).

Our results could appear to be at odds with the outcome of a study of two small low-metallicity (∼0.1 Z_⊙) star-forming regions in the extreme outer Galaxy (C. Yasui et al. 2009, 2010), where a small fraction of sources with NIR excess suggested very short disk lifetimes, even shorter than in the solar neighborhood. However, subsequent studies of other low-metallicity Galactic clusters indicate that the actual fraction of disk-bearing stars in the outer Galaxy is controversial (F. Cusano et al. 2011; C. Yasui et al. 2016) and subject to large age uncertainties (M. G. Guarcello et al. 2021). Nor is it yet clear what is the exact role of stellar feedback on the evolution and lifetime of circumstellar disks. M. C. Ramirez-Tannus et al. (2023) studied circumstellar disks in the massive Galactic cluster NGC 6357: while they were expecting the disks to be more affected by the UV radiation from massive stars than in the much smaller Lupus star-forming region, the authors found that the spectra in the two regions are quite similar.

Thus, from an observational point of view, there is currently no compelling reason to believe that at low metallicity circumstellar disks are short-lived. Furthermore, the prevailing theories requiring short lifetimes appear to be problematic for the formation of the planets known to exist around low-metallicity stars (S. Sigurdsson et al. 2003; A. Niedzielski et al. 2009; J. Setiawan et al. 2010), and might question the validity of the core-accretion model for the formation of planetesimals (J. L. Johnson & H. Li 2012). For planets to form, dust must settle in the midplane before the disk dissipates, but the settling timescale is longer at the low dust-to-gas ratios expected in these environments and it increases with distance from the star. If disk dissipation in low-metallicity stars were as quick as initially reported (C. Yasui et al. 2009), only very small rocky planets close to the star could form (J. L. Johnson & H. Li 2012).

In conclusion, our results indicate that, at the low metallicities typical of the early Universe, the disk lifetimes may be longer than what is observed in nearby star-forming regions, thus allowing more time for giant planets to form and grow than in higher-metallicity environments. This may have significant implications for our understanding of the formation of planetary systems in environments similar to those in place at Cosmic Noon.

We are very grateful to an anonymous referee for providing insightful and stimulating comments that have helped us improve the presentation of our results. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program JWST-GTO-1227. The specific observations analyzed in this work can be accessed via doi:10.17909/tatn-1114. M.M. and N.H. acknowledge support through a NASA/JWST grant 80NSSC22K0025, and M.M., N.H., and L.L. acknowledge support from the NSF through grant 2054178. O.N. acknowledges support from STScI Directorâs Discretionary Fund. K.F. acknowledges support through the ESA research fellowship. K.B. acknowledges support through Mini Grant INAF 2022 TRAME@JWST (TRacing the Accretion Metallicity rElationship with NIRSpec@JWST). O.C.J. acknowledges support from an STFC Webb fellowship. C.N. acknowledges the support of an STFC studentship (2645535). M.M., N.H., and L.L. acknowledge that a portion of their research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). A.S.H. is supported in part by an STScI Postdoctoral Fellowship. We thank Brunella Nisini, Teresa Giannini, Juan Alcalá, and Manuele Gangi for their discussions.

Facility: JWST - James Webb Space Telescope (NIRSpec, NIRCam).

Software: eMPT (N. Bonaventura et al. 2023), NIPS (C. Alves de Oliveira et al. 2018b).

We derived the accretion luminosities listed in Table 1 following the relationships of J. M. Alcalá et al. (2017) between accretion luminosity and line luminosity of the type $\mathrm{log}({L}_{\mathrm{acc}}/{L}_{\odot })=a\times \mathrm{log}({L}_{\mathrm{line}}/{L}_{\odot })+b$, where a and b are specific to each line and L_line is the line luminosity corrected for extinction. The values of L(Brγ) are those measured in this work, while those of L(Hα) are taken from G. De Marchi et al. (2011a). The coefficients (and their uncertainties) are a = 1.13(0.05) and b = 1.74(0.19) for Hα, and a = 1.19(0.10) and b = 4.02(0.51) for Brγ.

The extinction correction at the wavelength of Brγ was computed from the A_V values measured by G. De Marchi et al. (2011a) for each star (see Table 1) and converted to the K band by adopting the standard relationship (H. C. van de Hulst 1949; G. H. Rieke & M. J. Lebofsky 1985) A_K ≃ 0.1 A_V . The typical reddening is A_K = 0.05 ± 0.01 and after extinction correction L(Brγ) becomes about 5% brighter. Although the effect introduced by extinction is comparable to the typical calibration uncertainty, it is systematic and we corrected for it.

The mass accretion rates (Table 1) were computed from the freefall equation

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}\simeq \displaystyle \frac{G\,{M}_{* }\,{\dot{M}}_{\mathrm{acc}}}{{R}_{* }}\left(1-\displaystyle \frac{{R}_{* }}{{R}_{\mathrm{in}}}\right),\end{eqnarray} \tag{ A1 }$$

where G is the gravitational constant, M_* and R_* the mass and photospheric radius of the star, and R_in the inner radius of the accretion disk, from where the accretion flow is supposed to start. The values of M_* and R_* (listed in Table 1) are from G. De Marchi et al. (2011a). The value of R_in depends on how exactly the accretion disk is coupled with the magnetic field of the star and is therefore uncertain. Following E. Gullbring et al. (1998), we adopt R_in = 5 R_* for all PMS objects.

Concerning the uncertainties, besides those on the accretion luminosity, those on the R_*/M_* ratio also contribute to the uncertainty on ${\dot{M}}_{\mathrm{acc}}$. The uncertainty on R_*/M_* is dominated by the reddening correction (G. De Marchi et al. 2010). Given the typical A_V = 0.5 extinction of our sources, adopting a conservative 0.1 mag uncertainty on A_V results in an uncertainty of ∼3% on R_*/M_*. This is smaller than the systematic uncertainty on the mass alone, which G. De Marchi et al. (2010) estimate to be about 20% because the effects of the extinction on R_* and M_* are correlated. Therefore, the 0.25 dex uncertainty on the L(line) − L_acc conversion dominates the uncertainty on ${\dot{M}}_{\mathrm{acc}}$. The overall cumulative relative uncertainty on the mass accretion rate is estimated to be ∼0.42 dex (J. M. Alcalá et al. 2017).