© The Authors 2025

1 Introduction

The majority of stars in the Milky Way formed in clusters that contained at least one massive star (M_* ≥ 8 M_⊙) (Lada & Lada 2003). It is likely that our own solar system formed in such an environment (Pfalzner et al. 2015). Fundamental differences exist between the environments of massive and low-mass SFRs, in particular with respect to the intensity of the UV radiation field and the interstellar gas and dust density.

The majority of spectroscopy-based studies of PMS stars have looked towards low-mass SFRs, in part motivated by their proximity. These regions tend to be within a few hundred par-secs (pc), making it feasible to study the PMS stars with high signal to noise (S/N), and to spatially resolve features like outflows and protoplanetary disks. However, despite the advantages offered by these nearby regions, they do not represent the typical star and planet formation environment found in the Milky Way, or throughout the Universe. The lack of high mass stars, low stellar densities, and modest interstellar gas and dust densities are in sharp contrast to clustered star forming environments. The closest example of clustered star formation is in Orion, located at a distance of 414 pc (Menten et al. 2007). Orion, which contains several sights of massive star formation, most notably the Orion Nebula Cluster (ONC), is the most well studied massive SFR in the Milky Way (e.g. Ricci et al. 2008; Megeath et al. 2012; Lewis & Lada 2016; Berné et al. 2023). Even Orion represents a relatively modest cluster, with the ONC containing only a handful of stars with spectral type O, the most massive of which is θ¹ Ori C, an O6 star lying at the heart of the Trapezium cluster. In contrast, NGC 3603 contains upwards of 40 O-type stars, including several Wolf-Rayet stars, all located within an extremely compact volume in the centre of the cluster (Drissen et al. 1995). These massive stars inject tremendous radiative and mechanical feedback into their environment, ionising their surroundings to an extent that is several orders of magnitude more extreme compared to Orion, which is itself several orders of magnitude more extreme than low-mass SFRs. Like the ONC, there remain dense photodissociation regions (PDRs) consisting of ionised, atomic and molecular material surrounding the stars in NGC 3603. The impact that these prominent environmental factors have on star and planet formation remains largely unknown, though recently more attention has been given to the potential implications.

In O'dell et al. (1993), the authors reported on the discovery of a population of externally irradiated protoplanetary disks in the ONC, dubbed proplyds. The UV radiation from the massive stars was seen to be driving powerful thermal winds from the surface of nearby protoplanetary disks, with mass loss estimates of 10⁻⁶ M_⊙yr⁻¹. These mass loss rates would lead to disk dispersal in ∼10⁵ years, dramatically shorter than the dispersal timescales seen in low-mass SFRs ∼5 Myr (Mamajek 2009). In Haworth et al. (2018, 2023), the external photoevaporative mass loss rates of protoplanetary disks were calculated for a range of stellar masses, disk radii, and external UV field strengths. The UV field strength is given in Habing units G₀, with the typical interstellar medium having a value of 1 G₀. The authors found that in G₀ environments of ∼10⁴, comparable to the ONC and lower than in NGC 3603, disks around low and intermediate mass stars exhibit high mass loss rates down to disk radii of just 1AU. In Walsh et al. (2012), simulations were performed on an externally irradiated protoplanetary disk to understand the chemical changes that take place in the disk due to UV and X-ray radiation. They found that the UV radiation leads to a depletion of molecules in the inner disk (≤ 1 AU), and an enhancement of molecules outside of this radius. In Kuffmeier et al. (2023) the authors simulated star formation in regions of high gas and dust density, with column densities of interstellar material comparable to that found in NGC 3603 (Nürnberger et al. 2002; Schneider et al. 2015; Rocamora et al. 2024). The aim of this study was to investigate whether late-stage infall of interstellar material could rejuvenate protoplanetary disks, and possibly affect the evolution of the central star itself. They found that in such dense clustered environments not only does late stage infall commonly occur, but that in these episodes, on average >0.1 M_& is added to the star-disk system. For stars with a final mass of >1 M_⊙, more than 50% of their mass comes from late stage infall. In Winter et al. (2024), the authors reported a spatial correlation between the mass accretion rate and interstellar molecular gas density for a sample of young stars in Lupus, appearing to be in agreement with the results from Kuffmeier et al. (2023) that mass accretion can be influenced by the external environment.

This non-exhaustive list of studies highlights the need to consider the role of the external environment in the context of star and planet formation. Not only because there is compelling evidence that it can be a dominant influence on the evolution of planetary systems, but because these environmental influences are likely the norm, and hence affect the majority of stars and planets in the Galaxy.

With the launch of JWST, it is now possible to observe distant, massive SFRs in unprecedented detail. The sources observed in this study reside within the giant Galactic massive SFR of NGC 3603, at a distance of 7.2 ± 0.1 kpc (Drew et al. 2019). With its compact central cluster of OB stars, NGC 3603 is the Milky Way's closest analogue to a starburst cluster. In order to efficiently study the PMS stars in this cluster, we have utilized the Near Infrared Spectrograph's (NIRSpec) multi-object spectroscopy (MOS) mode. NIRSpec's MOS mode makes use of the novel Micro-Shutter Assembly (MSA), an array of over 250 000 microscopic doors that can be configured opened or closed, allowing for the simultaneous observation of up to hundreds of astrophysical spectra (Ferruit et al. 2022). We have obtained spectra of 100 sources in high resolution mode using the G235H – F170LP grating/filter combination.

In this paper, the spectra that we have classified as PMS stars are presented and discussed. In Section 2, we discuss the target selection. In Section 3, the data reduction process is explained. In Section 4, we describe additional calibration and processing steps including the nebular subtraction procedure that we have developed. In Section 5, we describe our methods including our spectral fitting procedure. In Section 6, we present our results including the spectral type, ages and masses and accretion rates of our sources. In Section 7, these results are discussed. In Section 8 our findings are summarised.

2 Selection of targets and observation strategy

The sources in this study were originally observed with Hubble Space Telescope's (HST) Wide Field Camera 3 (WFC3) (PID: 11360, PI: O'Connell, 2009). The filters employed in this observing programme were F555W, F625W, F814W, F110W and F160W. A different programme provided the F435W filter from HST/ACS observations (PID: 10602, PI: Maiz-Apellaniz, 2005). We have used these archival HST observations as well as our new NIRSpec observations to spectrally classify our sources (see Section 5.1). In Beccari et al. (2010), hereafter B10, the F555W and F814W filters were used to characterise the PMS population of NGC 3603. These sources were classified as either MS or PMS based on photometric H_a excess emission (with equivalent width EW > 10Å), which is known to trace accretion. Approximately 10 000 sources were observed in total. Of these 10 000 sources, 100 were selected to be observed by NIRSpec, of which 60 had been classified as PMS based on H_α excess above the 5 cc level. The remaining 40 did not show evidence of H_α excess above the 5 c level, and so were classified as MS stars. This strict selection criterion was employed in B10 to avoid the inclusion of sources with poorly subtracted nebular backgrounds or chro-mospheric emission masquerading as accretion, as both of these effects could lead to modest H_α excess emission. The approach of B10 was therefore highly conservative in identifying PMS stars. It is quite likely that many of the sources classified as MS are in fact PMS stars with intrinsically weak accretion signatures, or alternatively with significant NIR veiling, leading to a weaker detection of line emission. From our NIRSpec observations, we have classified 42 sources as PMS stars, based on the presence of at least two of the strong hydrogen emission lines Pa_α, Br_γ, and Brβ in emission above the chromospheric level, after performing corrections for photospheric absorption and veiling. Of the 42 sources that we have spectroscopically classified as PMS, 25 were classified as PMS in B10. The remaining 17 were classified as MS. Of the 17 sources in which our classification differs to B10, some of these sources display weak accretion signatures which could have been overlooked due to the strict selection scheme in B10. Additionally, a number of them lie in regions of complex nebulosity, requiring a careful subtraction of the nebular emission in order to properly reveal the underlying accretion related emission (see Section 4.2). However, there are a number of robustly accreting sources with bright emission lines in their NIRSpec spectra that had been classified as MS in B10. In these cases it is harder to reason why they may have been misclassified, though differences in methodology between this paper and B10 including background subtraction and EW measurement could play a role, along with intrinsic variability of PMS stars that strongly affects the strength of emission lines.

Compared to the sample used in our earlier work (Rogers et al. 2024c), the sources that we classify as PMS stars here include an additional 9 sources, bringing the total from 34 to 42. This is the result of improvements in the methodology allowing us to identify more sources with confidence, most importantly in the subtraction of the nebular background to reveal accretion related emission.

For our JWST NIRSpec observations, 3 micro-shutters were opened per source, forming a 'mini-slit'. The neighbouring shutters above and below the source shutter were opened in order to measure the nebular emission, and subtract it from the stellar spectra (see also De Marchi et al. 2024). This is a crucial procedure in a region like NGC 3603 where the nebular emission is bright. A 3 nod strategy was employed moving the source from the central shutter, to the upper shutter, and finally to the lower shutter. This pattern allowed for source spectra to be measured on different regions of the detector. The repeated exposures were eventually averaged together. This improved the S/N, while making it easier to remove cosmic rays, spurious pixels and other detector blemishes. The final observation strategy consisted of 3 unique MSA configurations, with 3 nods per configuration, resulting in 100 stellar spectra, and 600 nebular spectra.

In order to obtain a clean nebular spectrum for each target star, the micro-shutters observing the nebula must be free of stars. Due to the intense crowding of NGC 3603, this criterion placed tight constraints on which sources could be observed. The final list of targets was strongly influenced by how isolated each source was. Despite these precautions, it was not possible to completely avoid contamination of unwanted stars in nebular background micro-shutters. The number of usable nebular spectra decreased from 600 to 471 due to the presence of stars within the nebular micro-shutters. However, this was more than sufficient for the nebular subtraction.

3 Data reduction

3.1 NIPS

The uncalibrated data were largely reduced with the ESA Instrument Team's pipeline known as the NIRSpec Instrument Pipeline Software (NIPS) (Dorner et al. 2016). Additional reduction steps were also written specifically for these data which are discussed below. NIPS is a framework for spectral extraction of NIRSpec data from the count-rate maps, performing all major reduction steps from dark current and bias subtraction to flat fielding, wavelength and flux calibration, background subtraction and extraction, with the final product being the 1D extracted spectrum.

One of the final steps before extraction is the rectification of the spectrum, which is necessary because the dispersed NIRSpec spectra are curved along the detector. Rectification is performed in order to 'straighten' the spectrum. By doing this, it also resamples the spectrum onto a uniform wavelength grid (each wavelength bin Δλ is the same for every pixel). In our case, the rectified or "regular 2D" spectrum, was a count-rate map consisting of 3817 pixels in the dispersion direction and 7 pixels in the cross-dispersion direction.

The final data reduction step, namely the extraction, simply collapses the 2D spectrum by summing the 7 pixels along each column. Rather than working with the extracted 1D spectra, we chose to work with the rectified spectra, and performed further processing before eventual extraction in order to improve the S/N and aesthetic quality of the spectra.

3.2 Optimal extraction

During spectral extraction, a common approach is to simply sum over all of the pixels that the slit projects onto. In the case of a single micro-shutter, this is seven pixels along the detector column. For unresolved point sources, the majority of the flux is centred on just one or two detector rows. Using the simple sum extraction then leads to an unnecessary amount of noise being added to the spectrum, as most of the pixels contain little to no stellar light, but do contain detector noise. In order to minimise the noise that was added during the collapse of the spectrum, while including all pixels that contained stellar flux, an optimal extraction technique was employed, following the approach of Horne (1986). This method used all seven pixels per column to extract the spectra, but inversely weighted pixels based on their noise level. Including all pixels ensures that all of the stellar light is considered, while suppressing the pixels that contributed mostly noise.

A quick check was done on the S/N improvement achieved with optimal extraction. Two sources were chosen, one bright source with a V band magnitude of 16.2, and a faint source with V = 26.7. The two spectra were extracted using both the simple sum method, and the optimal extraction method. For the faint source, a region of the spectrum was defined that contained only continuum, and the S/N was measured to be 11.23. Using the optimal extraction technique, the signal-to-noise increased to 16.2, an improvement of 43%. For the bright source, the same approach yielded a signal-to-noise improvement of about 4%. This was not surprising, as optimal extraction was developed to suppress detector related noise. This is the dominant source of noise in faint spectra, and hence optimal extraction leads to a large improvement in those cases. Bright spectra are dominated by photon noise, and so optimal extraction has a less dramatic impact.

Finally, we combined the three nodded spectra for each source together. This was a multi-step process in order to again maximise the S/N while also removing spikes in the spectrum that had not been captured by the reduction pipeline. This procedure is described in detail in Appendix A.

4 Data corrections and calibrations

4.1 Absolute flux calibration

NIRSpec spectra are expected to have an absolute flux calibration accuracy of ∼10% (Gordon et al. 2022). We investigated whether this was consistent with our observations by comparing our PMS spectra with the HST photometry of our sources. These observations are all available on the Hubble Legacy Archive (HLA). We used the F160W (H band) flux of each source to assess the flux calibration. We performed an interpolation from the shortest wavelength of our spectra at 1.66 μm to the pivot wavelength of the F160W filter at 1.5369 μm. We did so by fitting either a 2D or 1D polynomial to the spectrum of each source. For the majority of our sources, we are observing the Rayleigh-Jeans tail of their spectrum, which is well approximated by a quadratic curve. For younger sources with flat or rising SEDs, the spectral shape around 1.66 μm is better approximated by a linear function, and so we used a 1D polynomial. We found that in order to bring the NIRSpec spectra into agreement with the F160W fluxes, a typical correction factor of just 1.02 needed to be applied. The 1σ variation of this scaling factor was 0.11. This is consistent with the 10% accuracy expected, indicating little to no variability in the continuum of our PMS source spectra between the HST observations taken in 2010 and the JWST observations in 2022.

4.2 Nebular background subtraction

4.2.1 The need for a scaled nebular subtraction

The subtraction of the nebular background from the source spectra proved to be a significant challenge. NIPS is capable of performing an automatic subtraction at detector level, in which it subtracts from the trace of the stellar spectrum, an equally sized trace containing nebular emission, as measured in the neighbouring micro-shutters. The performance of this subtraction routine was assessed in Rogers et al. (2022) based on simulated NIRSpec data for nebular emission with spatially uniform brightness. In that study, we found that the subtraction resulted in a residual nebular flux of 0.8% with a spread of 13%. After inspecting the real nebular spectra obtained with NIRSpec, we found that the nebular emission did in fact change in brightness by typically a few percent over the scale of a few micro-shutters (micro-shutter area = 0.46″ × 0.2″). The typical difference in emission line flux for Br₇ between neighbouring nebular regions was found to be 7%. In five out of 42 PMS sources, the corresponding nebular emission line fluxes differed by >20% between adjacent micro-shutters. In order to account for the spatially variable nebulosity, a scaling procedure was developed to subtract the nebular light from the stellar spectrum.

4.2.2 Can He I act as a scaling metric?

Following De Marchi et al. (2024), our brightness scaling procedure was based on the removal of the He I emission line doublet centred at 1.869am. This method assumes that the nebula is the dominant source of helium emission. Excitation from extreme-UV photons from the central OB stars is likely the dominant mechanism to produce this emission line, given the high excitation potential of this transition ∼23 eV. This doublet lies in a wavelength range that is very challenging to observe from the ground due to telluric absorption. It has also not been covered by the majority of infrared space telescopes, apart from the NICMOS instrument on HST, and now JWST NIRSpec. As such, this helium doublet has not been well studied or discussed in the literature. We have examined the available archival T Tauri and protostar spectra observed with NIRSpec, in order to determine whether this line may have a significant contribution from the circumstellar environment. The JWST programmes that observed these sources are given in Table 1, along with the name of the sources and the source type. All of these sources are located in Taurus-Auriga and Ophiuchus, where nebular emission is negligible or entirely absent. As such, the presence of He I 1.869 μm would strongly suggest an origin in the circum-stellar environment. We found that this He I line is absent in all protostar spectra, with an EW consistent with 0. In the case of the edge-on disks, we found weak He I emission in three of the five sources, with EWs ≤1 Å. In the remaining two edge-on disks, the feature was again absent. Additionally, this feature does not appear in absorption except for massive stars of spectral type mid-B and earlier (Husser et al. 2013). A number of the youngest and brightest sources in our sample, with the strongest emission lines, show no He I 1.869 μm emission in their spectra, even before nebular subtraction. These sources are so bright that the nebular recombination lines cannot be detected above the continuum level of the star. The lack of He I 1.869 μm in these sources demonstrates that even when hydrogen lines are strongly in emission, one should not necessarily expect to detect correspondingly strong circumstellar He I1.869 μm. At the same time, it seems more likely to find He I 1.869 μm in more active sources, as this line requires either very high temperatures or high degrees of irradiation. As such, it is unlikely to be produced by weakly accreting sources.

Our nebular subtraction routine works by scaling the subtraction until this He I line has been fully removed from the stellar spectrum. To quantify the potential over-subtraction caused by removing the He I line, we have re-run the nebular subtraction for five sources with high accretion rates, where the line could plausibly form in the circumstellar environment. To begin, we ran the nebular subtraction as before, removing the He I completely. We then measured the resulting line flux for Pa_α. We then performed the subtraction again, but stopped once the He IEW was reduced to EW = 0.5 Å. The resulting change in Pa_α flux was modest, increasing on average by 9% when the He I EW was reduced to 0.5 Å.

These results together suggests that the circumstellar environment for PMS stars can, at best, produce only weak He I 1. 869 μm emission, with EW at the level of <1 A. If this is the case for some of our more active sources, we may underestimate the accretion rate by of order 10%. More NIRSpec observations of CTTS and Herbig AeBe stars in low-mass star forming regions absent of nebular emission are needed to understand the extent to which He I 1.869 μm is produced in the circumstellar environment. For this study, we proceed by fully removing the line from our PMS spectra.

4.2.3 The scaling procedure

The procedure to scale our nebular subtraction worked as follows. The nebular spectrum was multiplied by a scaling factor, beginning with 0.01. The scaled nebular spectrum was then subtracted from the source spectrum. The He I EW in the source spectrum was measured. This was repeated with the scaling factor increasing in steps of 0.01 each time. The scaling factor that resulted in a He I EW as close to zero as possible was saved for that source, before moving on to the next source. In practice, after the scaled nebular subtraction, the stellar spectra had a typical He I EW of 0.011 ± 0.16 Å. The typical scaling factor was 0.94 ± 0.2. The scaling factor being slightly below unity with low scatter further supports that the He I 1.869 μm is nebular. Sources producing significant circumstellar He I 1.869 μm would require a larger scaling factor, which would be » 1. Outliers such as this are not seen in our sample.

As a result of the nodding pattern employed for our observations, we obtained six nebular spectra for each stellar spectrum (though some of these were not useable due to contamination). This meant that we could normally attempt the nebular subtraction six times per source, using a different nebular spectrum each time. This was done because we found that for a given source, some nebular spectra resulted in a cleaner subtraction than others. Here we define "clean" as being a subtraction in which there was minimal residual in the He I line, and none of the hydrogen lines became negative due to over-subtraction. The variation in performance of different nebular spectra was likely due to the poor pixel sampling. NIRSpec does not Nyquist sample the line spread function (LSF) at any wavelength for point sources using the R = 2700 mode. In many cases weak emission lines were sampled by just 1 or 2 pixels. This meant that even a small difference in the distribution of flux across those pixels from the nebular spectrum to the stellar spectrum could be enough to cause over-subtraction in one of the pixels, resulting in a negative flux. Higher sampling in combination with higher spectral resolution would spread the flux over more pixels, making the subtraction less sensitive to this. Having access to multiple nebular spectra per source allowed us to mostly circumvent this problem, and typically we could find a nebular spectrum that resulted in a clean subtraction. We have employed the conservative approach that if a clean subtraction was not possible for a source, it was classified as MS.

Fig. 1 Relationship between the EW of He 1 1.869 μm with: left, Paα; centre, Br7; and right, Br6.

4.2.4 Does removing He I remove H I?

Using the He I line as a metric for nebular subtraction relied on the assumption that the hydrogen recombination lines are also removed (i.e. there is a scaling relationship between nebular helium and hydrogen). In order to test the validity of this assumption, we measured the EW of the He I emission line, as well as the EW of the strongest H I emission lines: Pa_α, Br₇ and Br₆ for all of our nebular spectra. The resulting relationship is shown in Fig. 1, along with lines of best fit. All three hydrogen lines show a tight scaling relationship with He I. The dispersion around the line of best fit becomes larger for lines farther away in wavelength from He I. This implies that for longer wavelengths, other nebular lines may be needed to calibrate the nebular subtraction¹.

Equations (1), (2), and (3) show the coefficients of the best fitting lines for Pa_α, Br₇ and Br₆. The y-intercept of each equation indicates the EW of each hydrogen line when the EW of the He I line has been reduced to zero. In all three cases, some residual hydrogen emission is present according to these lines of best fit. A useful value to know is how much nebular flux may be left over in our subtracted stellar spectrum. In the case of Pa_a, there are 19.95 ± 2.26 Å of hydrogen emission remaining. These EW values are of course with respect to the nebular continuum, which is significantly lower than the stellar continuum. In order to convert this EW to a flux, we have multiplied it by the typical continuum level of our nebular spectra. This equated to a typical residual flux of 4.29 ± 0.49 × 10⁻¹⁷ erg/s/cm². For context, the median Pa_α flux from our PMS sources was 2.03 × 10⁻¹⁵ erg/s/cm². The residual nebular flux then corresponds to a typical under-subtraction at the level of ∼2%, which we do not deem significant with regard to the other uncertainties that affect our measurements. Given our argumentation in the previous section that we may even over-subtract some of our spectra, the true residual nebular flux is likely somewhere between +2% and −8% depending on the source in question. $$P{a}_{\alpha }=(15.36\pm 0.026)\times HeI-(19.95\pm 2.26)\AA$$(1) $$B{r}_7=(1.70\pm 0.003)\times HeI-(5.31\pm 0.31)\AA$$(2) $$B{r}_6=(5.04\pm 0.008)\times HeI-(-3.03\pm 0.77)\AA$$(3)

5 Methods

5.1 Determining spectral type, extinction and veiling

In order to establish the age, mass and accretion rate of our sources, we needed to first determine their spectral type, as well as their extinction, and excess continuum flux related to the protoplanetary disk and accretion. This excess flux is referred to as "veiling". Veiling is usually defined as the ratio of disk emission to photospheric emission, denoted by r_λ. We implemented a Markov chain Monte Carlo (MCMC) exploration procedure to simultaneously determine the stellar and disk parameters by

fitting our NIRSpec spectra and HST photometry with stellar models that had been both extinguished and veiled. Our method shares much in common with the approach outlined in Herczeg & Hillenbrand (2014) for optical spectra with moderate resolution. We also drew inspiration from Czekala et al. (2015) with the incorporation of a MCMC exploration to obtain meaningful uncertainty estimates on each parameter. Our approach consisted of first fitting the normalised NIRSpec spectrum of each source in order to obtain an estimate of T_eff and r_λ based on the minimum χ². These estimates would act as the basis of our priors for T_eff and r_λ. For the remaining parameters we used uninformative (flat) priors. Following this, we employed MCMC in order to determine more precise parameter values and their uncertainties.

5.1.1 Normalised spectrum fit

We opted to use the Phoenix stellar models from Husser et al. (2013) as our starting point. We degraded the spectral resolution of the Phoenix models to match NIRSpec using a Gaussian convolution and resampling onto a common wavelength axis. These models have three parameters, namely, effective temperature T_eff, surface gravity log(g) and metallicity [Fe/H]. We only considered a sub-sample of the Phoenix spectra, with 3000 < T_eff < 12 000, log(g) = 4.00 and [Fe/H] = 0.0. At the modest spectral resolution of NIRSpec, changes to absorption lines due to surface gravity are not significant, and so we elected to fix log(g) to the typical value for CTTS of 4.00 (e.g. Herczeg & Hillenbrand 2014). The metallicity of NGC 3603 has been investigated in other studies and found to be consistent with [Fe/H] = 0.00, and so we fixed it at this value. We review these studies briefly in Section 7.1. We normalised the NIRSpec spectra using a univariate spline function, setting the continuum level to unity. We ran a simple χ² fitting routine on each normalised spectrum, with 3000 < Teff < 12 000K and 0 < rA < 20. The values of the best fitting model, and the standard deviation of the ten best fitting models were used to generate the normally distributed priors for T_eff and r_λ during the MCMC part of this analysis.

5.1.2 Extinction

As our NIRSpec spectra cover a wavelength range that is not particularly sensitive to extinction, we also included HST photometry for each source. The filters employed were F435W (ACS), F555W, F625W, F814W, F110W, and F160W (WFC3). We performed aperture photometry for each source, accounting for the aperture correction and subtraction of the nebular background. This provided us with significantly better sensitivity to the extinction of each source. The optical fluxes are also sensitive to the photospheric continuum of the central stars, allowing for more accurate T_eff determinations.

In a previous work in which we examined the extinction towards NGC 3603 based on nebular emission line decrements (Rogers et al. 2024b), we used the NIR extinction curve parametrisation of Fitzpatrick & Massa (2009) (FM09). This extinction curve worked well with our nebular spectra, providing a parameterisation with a variable exponent a that produced better fits to our observations than other extinction curves with a constant exponent. However, the FM09 extinction curve does not cover some of our bluest optical data points, where we are most sensitive to extinction. Opting for an extinction curve with unbroken coverage from the optical to the infrared, we found that the recent extinction curve of Gordon et al. (2023) (G23) resulted in good fits to our observations when we tested it in our MCMC routine. The value of a used for the NIR portion of this extinction curve was α = 1.68467, which is consistent with the typical value that we determined in Rogers et al. (2024b), within the uncertainties.

In the future, if the FM09 curve is extended to include bluer optical wavelengths, it would be of great interest to compare it to G23.

We extinguished the Phoenix model spectra using the G23 extinction curve, setting R(V) = 4.8, based on our results in Rogers et al. (2024b). We did not allow this parameter to vary so that the number of free parameters was kept to a minimum. A(V) was permitted to vary between 0 and 10.

5.1.3 Excess veiling emission

Continuum emission from the protoplanetary disk can veil the stellar spectrum by reducing the apparent depth (i.e. the EW) of absorption lines. This emission is typically ascribed to hot dust at the dust sublimation radius of the protoplanetary disk (e.g. Millan-Gabet et al. 2006). One method to constrain NIR veiling in CTTS is to model the veiling emission as a black-body with a single temperature ranging from 500 to 2000 K (e.g. Muzerolle et al. 2003a; Cieza et al. 2005). This range represents the sublimation temperature for a variety of dust grain compositions and grain sizes expected in the protoplanetary disk. More recent works (e.g. Antoniucci et al. 2017; Alcalá et al. 2021) argued that hot gas in the inner disk also contributes towards veiling. In those studies, the veiling emission could only be successfully modelled with ablackbody spectrum with temperatures >2000 K. In McClure et al. (2013), three blackbody temperatures were used to fit the NIR veiling from 1 to 4.5 μm of CTTS: a cool blackbody ranging from 500 to 2000 K to represent the dust; a warm blackbody from 2000 K to the effective temperature of the star to represent optically thick gas in the inner disk; and finally the Rayleigh-Jeans tail of a hot blackbody of ∼8000 K to represent the accretion shock luminosity. The spectral wavelength range of McClure et al. (2013) is broader than our own, and more importantly goes to longer wavelengths where NIR veiling is stronger, justifying a multi-component fit to the veiling. For our observations we opted to fit a single blackbody, whose temperature could vary between 500 K up to the best fitting effective temperature of the star based on the normalised fit. This approximates excess emission from both the dust sublimation radius as well as gas in the inner disk without providing too much flexibility to overfit our observations. The value of r_A was again allowed to vary between 0 and 20.

5.1.4 Luminosity scaling term

In order to bring the Phoenix spectra into agreement with our observations, a luminosity scaling term was introduced. This was simply a multiplicative term that we applied across the entire spectrum. This term was varied until good agreement was found between the luminosity of the Phoenix spectrum and the observed spectrum. In general this term varied between 10⁻²³ and 10⁻²¹.

5.1.5 MCMC exploration

We employed a MCMC exploration with the python package emcee to determine realistic uncertainties on the best fitting parameter values. MCMC is an iterative sampling method to fit models to data. It works by generating a model based on the input parameters and randomly sampling values for each parameter from a pre-defined distribution (the priors). Each time a new model spectrum is generated, it is fitted to the observed spectrum, and the log-likelihood is calculated and taken as a goodness of fit metric.

If the current model provides a better fit than the previous model, then a new model is generated using the current parameter values as a starting point, plus a small random input. These small random steps in parameter space are accomplished by "walkers". If the current model is worse than before, then the next model is generated using the previous parameter values as the starting point, again, plus a small random input. This is repeated thousands of times until a reasonable range of parameter values has been determined, and no further exploration of the parameter space can improve the fit. This is known as convergence.

MCMC requires that arbitrarily small step sizes can be taken in parameter space when generating a new model to fit the observations. To enable this, we needed to interpolate the Phoenix model spectra. We used the radial basis function interpolation method from SciPy to do this. This method is computationally efficient, and works well for spectral features that change in a roughly linear fashion, which is a reasonable approximation for absorption lines with changing T_eff.

We allowed the MCMC algorithm to run for a maximum of 10000 iterations per walker, utilising 30 walkers in total. This meant that for each iteration, 30 new models were generated and fitted to the observed spectrum. This parallelisation allows for a more efficient exploration of parameter space. We checked for convergence by comparing the autocorrelation time of each parameter with the number of iterations that had been completed. If the number of iterations was >50 times the autocorrelation time, we stopped the procedure, as this indicates that the walkers had explored the parameter space sufficiently. Generally we reached convergence after approximately 3000-8000 iterations.

Figure 2 shows an example of a trace plot produced by the fitting procedure. This plot visualises how the walkers have explored the parameter space and on which values they have converged. Figure 3 is a corner plot, illustrating the distribution of parameter values. Figure 4 shows the source spectrum and the best fitting model spectrum. The remaining sources and their best fitting model spectra are available here.

Fig. 2 Trace plot of source 995.

5.2 Classifying continuum sources

A number of sources in our sample displayed extremely weak or no absorption lines in their NIRSpec spectra. Most of these sources exhibit rich emission line spectra, with transitions from H I, CO, Fe II among others. Strong emission line spectra are typically an indication of high levels of activity related to accretion. This serves to both fill in absorption lines with emission, while also leading to stronger NIR veiling, reducing the EW of absorption lines. Most of the continuum sources in our sample exhibit Class I or flat type SEDs in the NIR (see Section 6.2), which is also a sign of strong excess veiling emission from the protoplanetary disk, and suggests that the sources still retain substantial circumstellar material. Without photospheric absorption lines, we had to rely solely on fitting the SED of these sources in order to constrain the stellar parameters. The HST observations proved critical in this step. The stellar photosphere is well traced by optical photometry as the star's light dominates at wavelengths between 0.4 and 0.8 μm (Bell et al. 2013).

Fig. 3 Corner plot showing the distribution of parameter values.

5.2.1 Indirect clues about the spectral type of the continuum sources

The majority of the continuum sources are highly luminous and display no absorption lines in their NIRSpec spectra. The high luminosity of these sources can be explained if they are of intermediate to high mass, extremely young, or both. The complete lack of absorption lines from metals can be explained by one of three possibilities. The first, and least likely, is that NGC 3603 has a significantly sub-solar metallicity. This would render metal absorption lines intrinsically weak. However, numerous studies have examined the metallicity of NGC 3603, and point towards a solar, or near solar composition (e.g. Melnick et al. 1989; Lebouteiller et al. 2008). As such, we disregarded this possibility. Another option is that the metal absorption lines from the photosphere of our sources are so strongly veiled that they are no longer detectable. We have simulated how much veiling emission would be needed in order to render the strong metal absorption feature of Mg I at λ = 1.71 μm completely undetectable for a 4000K star, at the typical S/N of our continuum sources. We have assumed T_bb = 2000K. This would require a veiling factor of >20. At this level of veiling, the SED shape in the IR changes dramatically, as the blackbody spectrum of the veiling emission completely dominates over the stellar photosphere. We found that at this level of veiling, it was not possible to match the SED shape of our observations. As such, we do not favour extremely high veiling as an explanation. A Final explanation is that after ∼4000 K, metal absorption lines become weaker with increasing effective temperatures. At solar metallicity, metal lines vanish from spectra with temperatures of T_eff > 7000 K (Husser et al. 2013). These hotter photospheres are dominated by hydrogen absorption lines which, for our sources, are entirely filled in by emission lines, leaving the spectra with no detectable absorption. Veiling certainly plays a role in weakening the absorption lines further, but it seems more reasonable that a combination of moderate to high veiling with high effective temperatures is the reason why no metal absorption is seen for these five sources.

Fig. 4 Best fitting model spectrum shown in red, overplotted on top of the source spectrum shown in black.

5.2.2 Fitting the SEDs of our sources

We initially attempted to fit the SEDs of our sources with the Robitaille Young Stellar Object (YSO) SED models (Robitaille 2017; Richardson et al. 2024). Aside from the simplest family of models which represent a naked stellar photosphere, the Robitaille models attempt to fit numerous parameters related to the star itself, as well as the disk, envelope, geometry and inclination. We found that we lacked the wavelength coverage at longer wavelengths to constrain these model parameters. Given our modest wavelength range from 0.4 to 3 μm, we opted to only fit the optical HST photometry of our sources from 0.4 to 0.8 μm. By fitting only the optical observations, the underlying stellar properties could be determined without attempting to constrain the properties of the circumstellar environment. We performed a MCMC exploration, fitting the optical portion of our source spectra with the Phoenix stellar models, allowing T_eff, A(V) and the luminosity scaling term to vary. The models that we fitted are purely photospheric, while the observations show clear signs of photospheric emission plus disk emission. This necessitated that the photospheric models should not be significantly in excess of the observations at wavelengths >1 μm, where NIR veiling emission becomes strong.

We achieved good fits for our continuum sources. The best fitting temperatures tend to be relatively hot, with a median value of 6547 K compared to 4980 K for the rest of the sample.

5.2.3 Comparing SED fitting results to optical spectra

The central region of NGC 3603 was observed with MUSE by Kuncarayakti et al. (2016). A number of our sources fall within the field of view of those observations, providing us with a portion of their optical spectrum from 0.4 to 0.9 μm. Two of those sources, 354 and 152, are part of our continuum source sample, although they actually exhibit a single absorption line, Br_11-4. The presence and depth of this line indicates an effective temperature T_eff > 6000 K. The optical spectra of 354 and 152 are shown in Figures 5 and 6, respectively. The MUSE spectra overlap well with each source's HST photometry. Despite strong emission lines in their optical spectra, a number of hydrogen absorption lines are still detected, including Hβ, and a portion of the Paschen series. This provided us with the ability to compare the results of our SED fitting approach to a spectroscopic classification. We fitted the MUSE spectra of 354 and 152 with the Phoenix stellar models, allowing the temperature, extinction, and luminosity scaling term to vary, and calculating χ² each time. The resulting best fitting model for 152 has an effective temperature of 7200 K, compared to the SED best fitting temperature of 6185 K. The best fitting spectroscopic temperature for 354 is 7000 K, compared to 6875 K from SED fitting.

This shows that we underestimated T_eff in both cases, with the discrepancy being larger for source 152 than 354. What is clear from this exercise is that while our SED fitting results are more uncertain compared to a spectroscopic approach, we are more likely underestimating the true temperature of our sources, rather than overestimating it. This is consistent with the total lack of metal absorption lines in the continuum source spectra, which could only manifest for cool photospheres if an exceedingly high level of excess veiling emission were present. Our SED fitting shows that the continuum sources are relatively hot, intermediate mass sources, with a number of them belonging to the Herbig AeBe category. The stellar properties of these sources are shown in Table 2.

Fig. 5 Best fitting optical model spectrum for source 152 based on its MUSE spectrum. The black line corresponds to the MUSE source spectrum, the red line to the best fitting model spectrum, while the HST photometry is shown in grey.

Fig. 6 As above for source 354.

5.3 Measuring the recombination lines

To measure the recombination lines in the stellar spectra, we employed Monte Carlo simulations in order to propagate all of the uncertainties that arose from our processing steps and corrections. To do this, we generated 1000 realisations of each recombination line, allowing the lines to vary within their uncertainties for each realisation. We fitted a Gaussian profile to each realisation of the recombination lines using the non-linear least squares fitting routine curve fit from SciPy. We calculated the EW of the best fitting Gaussian, and converted the EW to a flux by multiplying it by the adjacent flux calibrated continuum. The final flux of each line is the median of the 1000 measurements. The uncertainty on the flux is the standard deviation of the 1000 measurements. The nebular recombination lines were measured in the same way. In the appendix, Figure B.1 shows the Pa_α emission line for all PMS sources, normalised to the continuum. The best fitting Gaussian profile is shown in each case as a red dashed line.

6 Results

In this section, we discuss the physical properties of the sources, derived from their spectra. The final values of the effective temperature, luminosity, extinction, and veiling were determined using our MCMC procedure as described above. The results for each source are given in Table 2.

6.1 Spectral properties

We have converted effective temperatures to spectral types using the conversion scheme of Pecaut & Mamajek (2013) for PMS stars. The majority of our sources (27/42) have spectral types G and F, with 5100 ≤ T_eff ≤ 6800 K. An additional 9 sources have type M and K with 3400 ≤ T_eff ≤ 5100 K. The remaining 6 sources have early F, A and late B types. The typical value of the extinction measured towards our sources is A(V) = 4.58 with a physical spread of ±1.85. This is consistent with the typical foreground extinction of A(V) = 4.5 found in previous studies (e.g. Melnick et al. 1989; Sung & Bessell 2004; Melena et al. 2008). The parameter that we were least able to constrain was the temperature of the NIR veiling blackbody emission. In many cases, the uncertainties are in the thousands of Kelvin. As we have discussed in Section 5.2, our lack of longer wavelength coverage means that our observations are not well suited to precisely determining NIR veiling parameters. We can in general conclude for a given source whether NIR veiling is present, but our ability to constrain the veiling properties beyond that is poor.

6.2 Spectral index

After correcting the sources for extinction, we measured the spectral index a of their NIRSpec spectra, defined as, $$\alpha =\frac{dlog\lambda F_{\lambda }}{dlog\lambda }$$(4)

from Wilking et al. (1989). This index is traditionally used to classify the evolutionary stage of PMS stars (i.e. Class I, II, III), with stars at earlier stages of evolution having a lower class. In Greene & Lada (1996), the authors defined as Class I sources those with α > 0.3, as flat spectrum sources those with 0.3 ≥ α ≥ −0.3, as Class II sources as those with −0.3 ≥ α ≥ −1.6, and as Class III sources those with α < −1.6. We classified our sources using the same scheme. The resulting histogram is shown in Figure 7. The majority of our sources (26/42) belong to Class III. We point out that, typically, the index a is measured in the range 1-10 pm using spectroscopy and/or photometry. Our wavelength range is therefore somewhat narrow, and we may underestimate the slope of our sources. CTTS SEDs in the IR are characterised by disk emission and dominated by hot dust longwards of 3 um, which results in a rising SED. Without access to these longer wavelengths where dust really "kicks in" (e.g. Furlan et al. 2011), we may be measuring a local minimum in the SED of our Class III sources. With this caveat in mind, it is interesting to find that the majority of our sources lack detectable continuum emission from their protoplanetary disk at 2-3 um (the H and K bands). Class III sources are thought to be essentially devoid of circumstellar material (Andre & Montmerle 1994), so our Class III sources would go undetected as disk-bearing PMS stars using only JHK broadband filter methods. However, the strong hydrogen emission lines in their spectra signal active and ongoing accretion, which necessitates that a gaseous disk still exists around the star. Accreting stars without detectable JHK excess have been observed in other Galactic massive star forming regions. In their study on the pre-main-sequence population of Trumpler 37 within the Cepheus OB2 association, Sicilia-Aguilar et al. (2005) reported that only 50% of their accreting stars showed excess emission in their JHK 2MASS photometry, despite showing strong H_α emission in their optical spectra. The authors go on to suggest external photo-evaporation as a possible reason for this. We discuss this as a possibility in Section 7.5.

6.3 Stellar ages, masses, and radii

In order to derive more information about the physical properties of our sources, we placed them in the Hertzsprung-Russel diagram (HRD), shown in Figure 8. We used the MIST stellar isochrones and evolutionary tracks of Choi et al. (2016) in order to determine their ages, masses and radii. The uncertainties on these parameters were determined via Monte Carlo error propagation by taking the 16^th and 84^th percentiles using a lognormal distribution to account for the fact that the age mass and radius should always be positive. The majority of the sources have ages between 0.1 to 3 Myr, typical for accreting CTTS. However, we also report a sub-sample of twelve stars consistent with ages > 10 Myr. Of those, four are consistent with ages > 15 Myr. Stellar masses and ages are listed in Table 2.

Fig. 7 Histogram of evolutionary classes for our sample.

6.4 Location of accreting sources within the cluster

We have plotted the location of the accreting PMS stars within NGC 3603 in Figure 9. Stars with ages of < 10Myr, taking into account their uncertainties, are shown in red. The remaining older stars with ages of >10 Myr are shown in blue. The older stars tend to be farther away from the cluster centre, while the younger stars exhibit more scatter, but are generally closer to the centre.

In order to test whether there is a significant difference between the projected distances of the old and young PMS stars, we have performed two statistical tests. First, we performed a Kolmogorov-Smirnoff (K-S) test, which measures the maximum difference (D) between the cumulative distribution functions of two samples in order to investigate whether they come from the same underlying population. The test produces two results, the K-S statistic itself D, and a p-value, which indicates whether D is the result of random chance. A p-value of <0.05 suggests that a significant difference exists between the two samples. We performed this test on the old and young stars by measuring their projected distances to the cluster centre. We obtained D = 0.550, p-value = 0.007. Although this result strongly suggests that the two samples are significantly different in terms of the distance to the cluster centre, our sample size is relatively small at just 42. In an attempt to improve the statistical robustness of our small sample size, we employed bootstrapping to generate confidence intervals for the median distance to the cluster centre for the old and young stars. Bootstrapping involves resampling the projected distances with replacement many times, calculating a median distance each time. This generates a distribution of median distances for both the old and young stars. The two distributions are then compared. The degree to which the distributions overlap indicates whether they are significantly different from Each other. The results from our bootstrapping analysis are shown in Figure 10. The younger stars have a median projected distance of 1.77pc, with a 95% confidence interval ranging from 1.49 to 2.04 pc. The older stars have a median projected distance of 2.57pc with a 95% confidence interval ranging from 2.11 to 2.80 pc. As these distributions do not overlap, this test along with the K-S test indicates that there is a meaningful difference in the projected distances of the old and young stars to the cluster centre.

6.5 Mass accretion rates

To compute the mass accretion rate Ṁ_acc of our sources, we first needed to calculate their accretion luminosity L_acc. To do this, we used the calibrations from Alcalá et al. (2017) (A17) for CTTS and Donehew & Brittain (2011) (D11) for Herbig AeBe stars to convert the line luminosity of Br₇ to L_acc. Having determined the masses and radii of our sources, we converted L_acc to Ṁ_acc using the formula, $${\dot{M}}_{acc}\sim 1.25\times \frac{L_{acc}{R}_{\ast }}{G{M}_{\ast }}$$(5)

from Gullbring et al. (1998); Hartmann et al. (1998). The emission line properties of our sources, including their Lacc and resulting MI_acc are given in Table C.1. The values of MI_acc span five orders of magnitude in our sample. Figure 11 shows Ṁ_acc of our sources with respect to M„. Unsurprisingly, we see that Ṁ_acc increases with increasing M„, though there is considerable scatter for the lowest and highest mass stars in our sample, well beyond their 1 cc uncertainties. There is no apparent trend between Ṁ_acc and M, for the low mass end of our sample. For M. > 1M_⊙, the accretion rate increases sharply with increasing stellar mass, reaching a peak of M_acc = 10⁻⁴ M_⊙ yr⁻¹. We fitted a line to the relationship between Ṁ_acc and M„, also shown in Figure 11. The best fitting coefficients to this line are $$log{\dot{M}}_{acc}=2.99(\pm 0.112)log{M}_{\odot }-8.0673(\pm 0.0525).$$(6)

We also include in Figure 11 the relationship found between Ṁ_acc and M. from D11. In that study, the authors measured the accretion rates of Herbig AeBe stars, but included a figure which combined observations from Herbig AeBe stars, intermediate mass T Tauri stars (Calvet et al. 2004), as well as CTTS (Muzerolle et al. 1998), spanning a range of masses from 0.112 M_⊙, corresponding to more than 3 orders of magnitude in mass. These measurements come predominantly from nearby low-mass SFRs. The relationship found for our sample of stars in NGC 3603 is significantly steeper than that of D11. This indicates that, for a given mass, the stars in our sample are accreting at higher rates than those from the combined sample taken from D11.

In Figure 12, we show the relationship between Ṁ_acc with age. As expected, MI_acc tends to decrease with age (e.g. Alcalá et al. 2014; Manara et al. 2015; Alcalá et al. 2017). The recovery of this relationship demonstrates that residual nebular emission does not significantly affect our measurements. We have fitted the relationship between Ṁ_acc and age. The line of best fit is shown in Figure 12. We have also included the line of best of fit reported in Sicilia-Aguilar et al. (2005), based on a combination of observations from several regions including Cepheus OB2, Taurus, Ophiuchus, Chamaeleon I, and TW Hydrae association (Muzerolle et al. 2000; Mannings et al. 2000). We have only included the portion of this line that is linear in log-space, starting at log age = 5.5. The Ṁ_acc versus age relationship found for our sample is offset and significantly less steep compared to the relationship reported in Sicilia-Aguilar et al. (2005).

We have demonstrated that our sample in NGC 3603 exhibits higher accretion rates for a given mass and age compared to nearby low-mass SFRs. We note that Cepheus OB2 is a massive SFR, harbouring the 06 star HD 206267. Sicilia-Aguilar et al. (2005) draw attention to several sources from Cepheus OB2 which also exhibit high Ṁ_acc and old ages. This behaviour is not seen in the observations from low-mass SFRs.

Fig. 8 HRD of the PMS sources using spectroscopically determined Teff and L*. The red points were classified spectroscopically. The orange diamonds were classified based on SED fitting. The lcr uncertainties are represented with red error bars. A selection of isochrones are shown as turquoise dashed lines. The ages are given in Myr on the right hand side.

Fig. 9 Location of accreting PMS stars in our sample. Sources with ages <10 Myr in red. Sources with ages >10 Myr in blue.

Fig. 10 Results of bootstrapping analysis. The projected distances of young stars in red. Old stars in blue. 95% confidence intervals shown as dashed red and blue lines. The median distance is shown as a dashed black line for each distribution.

6.6 Environmental correlations with M_acc

The mass accretion rates presented in Dil appear to follow a single relationship for the entire mass range 0.1-12 M©. Similar relationships have been found by Muzerolle et al. (2003b); Calvet et al. (2004); Herczeg & Hillenbrand (2008); Fang et al. (2009); Alcalá et al. (2014); Antoniucci et al. (2014); Manara et al. (2015), where a single power law can reasonably fit a wide range of stellar masses and accretion rates. There is significant scatter in the relationship that we have found between M_G and Af, for out sources. This motivated us to investigate whether there are environmental factors that correlate with the accretion rate, which may explain the scatter in Figure 11.

The environmental factors that we have considered are the ionised nebular gas emission, the molecular nebular gas emission, and the radial distance from the cluster centre. The ionised gas emission has been measured via the F656N H_α filter from HST WFC3. This was done by placing an annulus of r_in = 0.5″, r_out = 0″.6 around each source in the F656N image. These dimensions were chosen to ensure that light from the star had dropped off sufficiently with respect to the nebular emission. We calculated the median flux within this annulus, representing the typical nebular emission in H_α around each source. We have assumed that this emission is proportional to the density of the ionised nebular gas. The molecular gas emission was measured in the same way as H_α, using a H₂ narrow band image centred at 2.12 imi from the High Acuity Wide field K-band Imager (HAWK-I) at the Very Large Telescope. In this case we set r_in = 2″, r_out = 2.1″, owing to the larger PSF of HAWK-I. Again, we have assumed that the H₂ emission is proportional to the molecular gas density. The radial distances were obtained simply using the Euclidean distance formula, using the F814W HST image, which highlights the continuum emission from the stars.

In order to investigate whether these external factors may influence accretion, we first needed to remove the influence of both Af„ and age. As pointed out by De Marchi et al. (2011a, 2017), this is necessary because the mass accretion rate depends simultaneously on both parameters. The sources in our sample cover a wide range of ages and a PMS star of a given mass will have progressively smaller mass accretion rate as it ages. Following De Marchi et al. (2011a, 2017), we did this by performing a multivariate fit between Ṁ_acc, M_* and age t_*. The best fitting coefficients are, $$log{\dot{M}}_{acc}=1.682log{M}_{\ast }-0.795log{t}_{\ast }-\mathrm{2.731.}$$(7)

We then calculated the residuals of log Ṁ_acc and the best fitting line, and compared those residuals with the three environmental factors. We found no significant correlation between Ṁ_acc and ionised nebular gas emission, or the radial distance of the sources to the OB cluster. In the case of molecular gas emission however, we found a positive correlation with Ṁ_acc, as shown in Figure 13. The majority of the sample are associated with comparable, low levels of H₂ emission, finding themselves in the low density bubble near the central cluster. These sources exhibit a wide range of Ṁ_acc, indicating that a lack of high density molecular gas around a source does not imply a low accretion rate. However, in cases where accretion rates are higher than expected for a given age and mass, those same sources tend to be associated with regions of high density molecular gas, finding themselves in the prominent south-west and south-east pillars of the region.

Fig. 11 Accretion rate as a function of stellar mass. The size of each circle is proportional to age. The largest circle corresponds to 19 Myrs. The best fitting line to our observations is shown in brown, while the best fitting line from Dll is shown in black.

Fig. 12 Accretion rate as a function of stellar age. The size of each circle is proportional to M⊙ The largest circle corresponds to 7 MQ. The best fitting line to our observations is shown in brown. The best fitting line from Sicilia-Aguilar et al. (2005) is shown in black.

Fig. 13 Normalised accretion rate as a function of the nebular H2 brightness around each source. Sources associated with above average levels of nebular H2 brightness are coloured red.

6.7 Revisiting the $L_{P{a}_{\alpha }}$ relationship

In Rogers et al. (2024c), we established the first relationship between L_acc and $L_{P{a}_{\alpha }}$. Since the completion of that work, we improved our methodology related to correcting for extinction, spectral classification, flux calibration and nebular subtraction. We present here an updated relationship. The net result of our changes and improvements in methodology is a significant reduction in the uncertainties of the best fitting line between L_acc and $L_{P{a}_{\alpha }}$, with no changes to the best fit coefficients. $$log{L}_{acc}=1.42(0.08)log{L}_{P{a}_a}+3.33(0.17).$$(8)

In Figure 14, we show the relationship between the line fluxes for Br₇ and Pa_α. Owing to our high quality observations and improved methodology we were able to measure both of these lines with high precision, leading to small statistical uncertainties for each flux measurement. To calculate L_acc, we converted our Br₇ line fluxes to luminosities, and then used the calibrations of A17 and Dll to convert to L_acc. The coefficients of the revised relationship are identical to Rogers et al. (2024c), but with the uncertainties reduced by a factor ∼2.5. Figure 15 shows the resulting relationship between Pa_a and L_acc. The best fitting line coefficients are given in Equation (8). What is clear by comparing Figures 14 and 15 side by side is that the dominant uncertainty in this relationship is inherited from the conversion to L_acc from A17 and Dll. This is a result of the more challenging near-UV shock modelling approach that is used to directly determine L_acc in A17 and Dll. This approach unambiguously probes the material being accreted onto the star, but high extinction in the UV as well as absorption by Earth's atmosphere necessarily leads to larger statistical uncertainties.

Fig. 14 Relationship between the line flux of Br7 and Paα. The brown line represents the line of best fit and the grey shaded area shows the uncertainty of the fit.

7 Discussion

In this section, we discuss and interpret our results from above. This includes a brief review of the metallicity of NGC 3603 from previous studies, a discussion of disk tracing features detected in the NIRSpec spectra as well as the conspicuous lack of circum-stellar H₂ emission. The impact of external photoevaporation on our sources is considered and the ages of our sources are discussed and compared to stars in low-mass SFRs. The possible implications of the environmental correlation between Ṁ_acc and nebular H₂ emission are considered.

7.1 The metallicity of NGC 3603

We opted to fit our source spectra with only solar metallicity models corresponding to [Fe/H] = 0.0. The metallicity of NGC 3603 has been investigated before, typically by obtaining nebular spectra of the H II region, and using recombination lines from metals like Ne and O to infer a metallicity (e.g. Melnick et al. 1989; Lebouteiller et al. 2008). These studies report that NGC 3603 appears to be consistent with solar metallicity. As such, we have a strong prior belief that the metallicity of the stars should also be consistent with solar metallicity. PMS spectra are not well suited to a rigorous metallicity investigation, despite there being a number of strong metal absorption lines in our wavelength range from Al, Mg, Ca and Na. This is because the strength of metal absorption lines can be reduced due to NIR veiling from the protoplanetary disk. If the veiling is not fully characterised, it can lead to the appearance of a low metallicity due to apparently weak metal absorption. The metal species mentioned are most strongly in absorption for sources with low effective temperatures 3500 <T_eff< 5000 K. At higher temperatures, these metals can no longer survive in the stellar photosphere. For mid-G type stars and earlier at solar metallicity, metal lines are typically no longer visible (Husser et al. 2013). In the case of MS stars, it is therefore only the stellar temperature and metallicity that can affect the metal line absorption depths. Conversely, in the case of PMS stars, an additional variable is veiling, which itself is challenging to fully characterise. A possible solution to remove the influence of veiling in determining the metallicity is to look for metal lines, close to each other in wavelength, whose ratio changes with metallicity. As the lines are near each other, they experience approximately the same amount of veiling, and so their ratio should not change due to veiling. We did find one such ratio for the Al lines at 1.67189 and 1.67505 μm. Their ratio inverts when going from [Fe/H] = 0.0 to [Fe/H] = +1.0. However, none of our sources were well fitted with super-solar metallicity, and so this ratio was not used. We did attempt to fit a number of sources with sub-solar metallicity models, at [Fe/H] = −1.0, but the goodness of fit was comparable to the solar metallicity model best fit. Therefore, we have no reason to believe that the metallicity of our sources in NGC 3603 is different from the solar value.

Fig. 15 Relationship between Lacc and Paα. The brown line represents the line of best fit and the grey shaded area shows the uncertainty of the tit.

7.2 Hydrogen recombination lines

All 42 PMS spectra display some or all of Pa_α, Br₆ and Br₇ in emission. Pa_α A 1.875 μm and Br₆ λ 2.662 μm are two strong NIR lines that fall within wavelength ranges that are not typically recoverable from the ground due to telluric absorption. While it is possible to attempt to remove these telluric features with post-processing (e.g. Vacca et al. 2003; Kausch et al. 2015), the wavelength ranges of 1.8-1.9 μm and 2.5-2.7 μm are often omitted in NIR spectroscopic studies. JWST NIRSpec provides us with access to both of these lines simultaneously for the first time. This has facilitated the calibration of both of these lines as accretion tracers Rogers et al. (2024c). A subset of sources display additional Brackett lines as well as Pfund lines in emission. These sources are more massive and younger than the sample average. Both of these characteristics are typically accompanied by higher accretion and or ejection activity, which likely explains why these sources feature more line emission. The Pfund series is ∼3 times weaker than the Brackett series, making it observa-tionally challenging to detect with high S/N. Pfund lines have been detected with ground based observatories in star forming regions within <1000 pc (Maaskant et al. 2011; Koutoulaki et al. 2018), as well as in more distant HII regions from massive young stars (Bik et al. 2006), and have been detected in CTTS stars in the Small Magellanic Cloud Jones et al. (2022). In these cases only a limited selection of the lines are visible, again due to telluric absorption. The uninterrupted Pfund series has been studied with space based telescopes such as the Infrared Space Observatory (ISO) Benedettini et al. (1998). Due to the relative weakness of the series as well telluric absorption, the Pfund series is not as well studied as other NIR or optical hydrogen series. In Salyk et al. (2013), the authors showed that Pf₈ (λ 4.6538 μm) can be used to measure the accretion rate of PMS stars, making it a useful diagnostic for deeply embedded objects as extinction at this wavelength is extremely low compared to optical wavelengths. In Rogers et al. (2024a), the full width at half maximum (FWHM) of Pfund and Brackett emission lines was used to kinematically infer that these lines have an origin in the magnetospheric accretion flow, rather than from winds for a sample of intermediate mass stars, including a number of Herbig Ae stars. The Pfund series is shown in Figure 16 for a number of sources.

7.3 Disk tracing species

7.3.1 CO bandheads

Disk tracing species Fe II 1.688 μm and CO bandheads are detected in emission in a number of our PMS spectra. CO bandhead emission is thought to arise from within the dust sublimation radius of the protoplanetary disk, at relatively high temperatures (2000-5000 K) and densities (n_H > 10¹⁰ cm⁻³). An origin in the protoplanetary disk has been supported by studies at high spectral (Ilee et al. 2013; Bik & Thi 2004) and spatial (O'Garatti et al. 2020; Koutoulaki et al. 2021) resolution. It is thought that self-shielding occurs at these high densities, which protects the bulk of the CO from dissociation via UV photons from the central star (and in the case of our sources in NGC 3603, from the OB cluster). The profile of the CO bandheads is typically composed of a 'blue shoulder' due to doppler broadening of the line (Ilee et al. 2013). The red wing of the line is composed of a series of closely spaced vibrational transitions, known as J-lines. The J-lines are not resolved at NIRSpec's spectral resolution, giving the bandhead the appearance of a long red tail.

Of the 42 sources in this study, 5 show CO bandhead emission (12%). This is below the detection rate of ∼25% found in Ilee et al. (2013). These sources are 152, 185, 238, 251, and 354. All of them are of intermediate mass, ranging from 4 to 7 M_e. For sources 152 and 354, the first 6 bandheads are clearly detected. For both of these sources, there is tentative evidence of the blue-shoulder, especially in the strongest transition CO2-0-In the cases of sources 185, 238 and 251, the bandheads are noticeably weaker in terms of their EW, and no blue shoulder is detected. The weaker appearance of their bandheads is likely a result the stronger stellar continuum, as these sources are twice as massive as 152 and 354. This leads to veiling of the band-heads, reducing their apparent strength. All sources with CO bandhead emission also display strong hydrogen recombination lines and exhibit significant NIR excess, to the extent that no absorption lines are detected in their spectra. Figure 17 shows the rectified CO bandhead spectrum of each source, with an offset on the y-axis for clarity.

Fig. 16 Five PMS sources that display Pfund recombination lines. Their spectra have been normalised to unity, and offset in this figure for visual clarity. The Pfund lines are marked with a green dashed line. Br6 is also present, marked with a blue dashed line. Note that a number of pseudo absorption features are artificially created as a result of our normalisation.

7.3.2 Fe II

The Fe II fluorescent emission line at 1.688 μm is detected in just three of our sources: 251, 469, and 823. This line is more commonly detected in intermediate and high mass PMS stars (Porter et al. 1998), as it is likely pumped by Lyman photons, and hence a strong UV source is required. This is true for our sample of stars, with the line being detected in some of our most massive sources. As Fe is easily ionised to high ionisation states, in order for Fe II to survive, it must reside in dense gas. It is thought to trace dense material of the inner disk, in a region outside of where the CO bandheads are produced (Lumsden et al. 2012), and is known to originate in the inner disk of classical Be stars (Carciofi & Bjorkman 2006). As such, its origin in the protoplanetary disk of PMS stars seems reasonable. This line is located directly beside Br₁₁, and is notably narrower than Br₁₁ in all three sources. This supports the idea that Fe II is not produced in an outflow or accretion flow, as the high velocities involved would cause considerable line broadening. The presence of CO bandheads and Fe II in conjunction with strong continuum NIR excess indicates that sources with sufficient mass can retain a significant gas rich inner disk while being strongly irradiated by massive stars for several 10^s years.

Fig. 17 The normalised spectra of the five CO bandhead emission sources. The spectra have been offset for visual clarity.

7.4 Lack of circumstellar H₂

Despite our large, high S/N sample of PMS sources of varying masses and stages of evolution, we did not detect any H₂ emission from any of our sources. The following discussion is largely qualitative in nature, and is not intended as a robust explanation to the lack of circumstellar H₂, but rather to consider some potentially important physical mechanisms that could be at play for our sources.

H₂ is commonly seen in emission from PMS stars, coming either from the protoplanetary disk itself (Weintraub et al. 2000; Bary et al. 2002) or, more commonly, tracing a wind or outflow (Takami et al. 2007; Beck et al. 2008; Harsono et al. 2023). We detect numerous H₂ emission lines in the nebular spectra. The nebular spectra that show these lines most prominently are associated with the dense pillars of material seen in the south-west and south-east of NGC 3603. However, the absence of circumstellar H₂ could be due to a number of reasons.

Given the large distance to NGC 3603, emission from the protoplanetary disk or possible outflows is spatially unresolved out to distances of ∼700 AU. Emission from the star, disk and possible outflow are therefore all contained within the same beam. If the emission from H₂ is intrinsically faint, it may be too weak to detect above the continuum level of the stars. On the other hand, the robust detection of CO bandheads, as well as other disk-tracing species such as Fe II 1.688 μm likely rules this explanation out. When these features are detected in the spectra of young stars, H₂ tends to also be present, and is generally as strong or stronger than the CO or Fe II (e.g. Davis et al. 2011).

Another possible explanation is that, if the bulk of the H₂ emission would come from a wind or jet as is typically found for CTTS and Herbig AeBe stars, this outflow would experience external photoevaporation from the central OB stars. This could dissociate or even ionise the outflows. Externally ionised Herbig-Haro jets have been observed from numerous young sources in Orion (Reipurth et al. 1998; Bally & Reipurth 2001, 2003). Those observations spatially resolved the jet from the central star, and unlike typical molecular Herbig-Haro outflows, the spectra of those jets resemble nebular emission spectra from photoionised gas. It is therefore likely that a similar, or possibly even enhanced form of this jet ionisation occurs also in NGC 3603, given its much larger population of early O type and Wolf-Rayet stars compared to Orion. In Reipurth et al. (1998), the authors even suggest that launching of the jet may be entirely quenched on the side that faces the OB cluster, as the four sources in their study show significantly less developed jet lobes on the directly irradiated side, while the more prominent lobes face away, possibly shadowed by the protoplanetary disk itself.

The recent work of Kirwan et al. (2023) partially answers this question some twenty years later. Here, the authors studied the enigmatic Orion proplyd 2AA-AA0, which launches a prominent bipolar jet that is apparently unaffected by the UV environment. They reason that the jet was likely launched after the onset of external photoevaporation, and is shielded from much of the UV radiation by the envelope of the proplyd itself. This protection can only last until the jet reaches the ionisation front and penetrates the cusp of the proplyd envelope, after which it is exposed to the external environment. It seems plausible then that in NGC 3603 any H₂ outflows could have already been destroyed given the four orders of magnitude more intense UV field compared to Orion (Rollig et al. 2011), which would result in an ionisation front closer to the disk surface Winter & Haworth (2022).

This is consistent with the work of Reiter & Smith (2013) and Reiter et al. (2016), who used HST narrow band imaging of H_a and [Fe II] 1.644 μm (not to be confused with the permitted Fe II 1.688 μm) to detect and characterise externally ionised outflows from young stars in the Carina star forming association near the massive SFR Trumpler 14. Trumpler 14 contains at least 74 O stars (Berlanas et al. 2023). The existence of prominent, highly collimated jets in this region demonstrates that jet launching need not be entirely switched off in the presence of high UV radiation. These jets are consistent with the externally ionised jets found by Reipurth et al. (1998), producing predominantly atomic and ionic emission lines. In a more recent work using the IWST Early Release Observations (ERO) of NGC 3324, adjacent to the Carina nebula, Reiter et al. (2022) detect outflows traced by both Pa_α as well as H₂, with the vast majority of the H₂ outflows being found within molecular clouds where they are presumably protected from the harsh UV environment.

In summary, we suggest that the lack of circumstellar H₂ detected in our sample may result from external ionisation of any initially molecular outflows by the central OB cluster, as has been seen in other massive star forming regions.

7.5 External photoevaporation in NGC 3603

External photoevaporation is the defining process associated with the objects known as proplyds. Proplyds are characterised by their unique teardrop shaped cocoon, enveloping the protoplanetary disk and star, formed from material boiled off of the disk surface. For NGC 3603, we lack the spatial resolution to classify our sources as proplyds based on this criterium (although based on the total lack of circumstellar H₂ emission, it is highly unlikely that our sources are similar to the proplyds found throughout Orion). Regardless, our sources can unambiguously be classified as externally irradiated stars. As such, external photoevaporation has likely played a role in the evolution of these sources at some point in their lifetime.

External photoevaporation is a possible explanation for the low levels of NIR excess that we have measured towards the majority of our sources. This process truncates the protoplanetary disk, eroding the disk from the outside in (e.g. Ndugu et al. 2024). A smaller disk results in a smaller emitting area, and hence less NIR excess emission. Disk sizes in the Orion Nebular Cluster have been shown to be smaller than in other SFRs and even to scale with the distance to the central ionising OB stars (Boyden & Eisner 2020). With that in mind, most of the NIR excess is thought to originate from the inner disk within ∼ 1 AU of the star. Material this close to the star is located deep in the potential well, and is therefore more difficult to pho-toevaporate. The disks in NGC 3603 are more strongly irradiated compared to the disks in Orion. They will therefore have higher sound speeds due to higher temperatures. This serves to reduce the gravitational radius (Winter & Haworth 2022), which could allow for inner disk depletion. It has also been shown in Haworth et al. (2018) that in high G₀ environments, even the inner disk can be significantly impacted by external photoevaporation.

Without exact knowledge of the UV field at the surface of the disk, it is unclear whether a significant fraction of the inner disk could be lost due to external photoevaporation, and whether this could explain the lack of excess emission at 2-3 um. Several of the results that we have presented here also seem to contradict significantly depleted disks. The disk lifetime is expected to be dramatically reduced due to external photoevaporation (Henney & O'dell 1999), but we have determined that many of our sources are still accreting after 10 or even 20 Myr. Additionally, the accretion rates of our sources are higher than stars with comparable masses and ages based on other studies in the literature. This has already been been pointed out by De Marchi et al. (2017) from the Ha excess of sources in other massive star forming regions. Both of these results suggest that a robust gaseous disk survives around our sources. As such, it remains unclear how external photoevaporation has impacted our sources. Observations at longer wavelengths, for example with JWST's MIRI, would directly probe more complex chemical species in the disks. These molecular features are expected to change significantly as a result of UV driven photochemistry (Walsh et al. 2012, 2013), though recent results from MIRI seem to contradict this Ramfrez-Tannus et al. (2023). Follow up observations at longer wavelengths could also reveal to what extent the disks may have been truncated due to external photoevaporation.

7.6 How long can accretion last for stars in NGC 3603?

In low-mass SFRs with low levels of UV radiation, disk dispersal timescales are typically around 5 to 10 Myr (e.g. Armitage et al. 2003). Long lived disks with ages of 40 to 50 Myr have been found around very low mass stars, so called "Peter-Pan" stars. This has been attributed to below average levels of UV radiation in their environment (Coleman & Haworth 2020), resulting in lower than average levels of photoevaporation. It is therefore surprising that protoplanetary disks undergoing active accretion could survive for ∼15 Myr in an environment like NGC 3603, which has some of the highest levels of UV radiation in the Galaxy (Melena et al. 2008). The disk fraction in such an environment is expected to be close to 0% after 7 Myr (Guarcello et al. 2010).

This result of long lived disks in NGC 3603 is in agreement with the photometric study of B10, who found evidence of an old PMS population still undergoing accretion based on H_a excess emission using HST narrow band imaging. In that study, a single average extinction correction was applied to all sources of A(V) = 4.5, using the extinction curve of Cardelli et al. (1989) with R(V) = 3.55. We have determined the extinction spectroscopically for each source individually, and have used a higher value of R(V) = 4.8. The average value of A(V) that we found was 4.04, with a dispersion of 2.1 across the sources. The lower value of R(V) and higher value of A(V) used in B10 tends to push the sources towards bluer colours, making them appear older and closer to the zero-age-main-sequence (ZAMS) line. Indeed, if we apply the extinction correction to our sources as was done in B10, our old sources move from ages of ∼15 Myr to ∼50 Myr. Despite these differences in methodology, our spectroscopic findings are in general agreement with B10, that accreting PMS stars survive in NGC 3603 well beyond the expected cut-off point of ∼7 Myr.

Disk dispersal timescales have been inferred observationally by searching for stars with NIR excess emission in star forming regions with different ages (Strom et al. 1989; Haisch et al. 2001; Hillenbrand 2005; Yasui et al. 2009). Searching for NIR excess could be done photometrically with J, H, K bands, making it an efficient method to probe disk lifetimes in many different star forming regions. We have found spectroscopically that 26 of our 42 PMS stars show no detectable NIR excess in the J, H, K bands, but do show strong hydrogen recombination lines indicating ongoing accretion. We argue that emission lines provide a more sensitive means of inferring ongoing accretion, as these lines more directly probe accretion related activity compared to NIR excess, and are easier to detect in spectra compared to weak excess emission associated with the later stages of star formation.

Thanathibodee et al. (2018) studied three sources with low mass accretion rates from the UV to NIR. Two of the sources showed unambiguous evidence of ongoing, albeit weak, accretion based not only on hydrogen emission lines, but also red-shifted absorption in the NIR He I 21.083 um line, tracing infalling material. These sources showed negligible NIR excess, and a modest reservoir of remaining gas in the inner disk based on marginal detections of FUV H2 emission. More sensitive facilities like JWST can detect faint spectroscopic signatures of accretion that were not captured by earlier efforts using broadband NIR excess.

A recent paper from Winter et al. (2024) (W24) has examined whether PMS stars in regions of enhanced interstellar gas and dust density could experience disk replenishment by infall from the interstellar medium. They found a statistically significant positive correlation between the mass accretion rate of the sources, and the gas density around those sources. This result demonstrated that a spatial correlation exists with the accretion rate of stars. Infalling streams of cold molecular material have been detected towards numerous PMS stars (e.g. Tobin & Sheehan 2024). Their detection suggests that disk replenishment is a feasible mechanism to extend the lifetime of protoplanetary disks. Late stage infall could help explain why protoplanetary disks in massive SFRs can survive for longer than expected, even under the influence of external photoevaporation.

In a recent study, De Marchi et al. (2024) presented a sample of old (∼30 Myr) accreting PMS stars in the massive SFR NGC 346 in the Small Magellanic Cloud using NIRSpec MSA observations. That study employed a similar approach to determining the age of their stars as we have done in this study. Those authors have studied the Magellanic Clouds extensively via HST photometry (Sabbi et al. 2006; De Marchi et al. 2011a, 2013, 2017; Biazzo et al. 2019; Tsilia et al. 2023), and have consistently found populations of accreting stars with ages well beyond what is thought possible based on current disk evolution theory. We note for the reader that many of these same authors conducted the survey presented in B10. The observations that we have presented here as well as those presented in De Marchi et al. (2024) serve as spectroscopic follow ups to those photometric studies that claimed to have found old, strongly accreting PMS stars based purely on photometric measurements. These follow up observations prove that the results inferred from photometry were largely accurate, that disks around PMS stars can survive for much longer than we think, and indicate that there is are significant gaps in our understanding about how stars and protoplanetary disks form and evolve in massive SFRs.

Combining the high sensitivity of NIRSpec with the use of emission lines as disk tracers has enabled the detection of accreting PMS stars that would have largely gone undetected based solely on NIR continuum excess. Disk fractions inferred purely from NIR continuum excess are likely underestimates. We have shown that most of our sources show no detectable excess in the J, H, K bands, but still have high accretion rates.

7.7 How accurate are stellar isochrones?

The use of stellar isochrones and evolutionary tracks are commonly employed in star formation studies to determine stellar properties like age and mass. Mass determinations from PMS evolutionary tracks are relatively reliable, with accuracies of ∼10% (Stassun et al. 2014). Age estimates on the other hand typically carry larger uncertainties of a factor 2-3 (Soderblom 2010). In Bell et al. (2013), the authors re-assessed the ages of 13 star forming regions by fitting the MS population with stellar evolutionary models, and fitting the PMS population with their own semi-empirical PMS isochrones, finding good agreement between the two methods, which each rely on different stellar physics. The ages determined by the authors were typically ∼2 times older than previous estimates in the literature, indicating that stellar isochrones tend to underestimate the true age of PMS stars. If this can be applied to the MIST isochrones that we have employed, then our stars with ages of ∼15 Myr will in fact be closer to ∼30 Myr in age. Indeed, the point is raised by Bell et al. (2013) that the timescale of star formation may be longer than we think due to a systematic error in isochronal ages. This systematic error serves to push stars towards younger ages. This doesn't change the fact that the timescale for disk dispersal appears relatively consistent across numerous different isochronal models for many star formation regions, even if this agreement is only consistent on a relative level. Studies have found that after ∼5 Myr, only 20% of stars are still surrounded by a protoplanetary disk (e.g. Carpenter et al. 2006; Dahm & Hillenbrand 2007), with this fraction being apparently even lower in massive SFRs (Guarcello et al. 2021). As pointed out by De Marchi et al. (2024), stars with active accretion signatures in Galactic low-mass SFRs are not found in the region of the HRD where our oldest PMS stars are located.

While the ages of stars as determined by isochronal fitting may be systematically offset and inaccurate on an absolute level, the relative difference in ages demonstrates that the timescale of protoplanetary disk dispersal differs between stars in low-mass SFRs and massive SFRs.

7.8 Systematic uncertainties affecting the HRD

Systematic uncertainties also affect the position of stars in the HRD, such as uncertainty in the distance to the region as well as the intrinsic spread of distances to sources within the cluster. In a review of Scorpius-Centaurus, the nearest OB association (Preibisch & Mamajek 2008), the authors point out that the physical spread of distances in Upper Scorpius is ±20 pc, which is a relatively large fraction of the total distance to the region at 145 pc. This introduces a significant degree of uncertainty to the positions of the sources on the HRD. In regards to NGC 3603, this is less significant for two reasons. First, the distance to NGC 3603 is large, meaning that the intra-cluster distances would also need to be large in order for it to become an important source of uncertainty. Second, while Upper-Scorpius is a large association, NGC 3603 is a dense, compact cluster, with a cluster radius of ∼5 pc (Nürnberger & Petr-Gotzens 2002). Based on this, the small differences in distance to the individual sources in our sample should not significantly affect their position on the HRD. A consideration that has not been discussed so far is multiplicity. It is likely that more than half of our sample are in binaries or higher order systems (Reipurth & Zinnecker 1993). The presence of companions serves to increase the total luminosity of our sources. By not taking this into account in our analysis, we have undoubtedly introduced additional uncertainty into their positions on the HRD. However, given that this effect serves to push our sources up, away from the ZAMS, the age estimates are likely underestimated, rather than overestimated.

7.9 Spatial distribution of PMS stars in NGC 3603

A significant difference at the 2o- level in the projected distances to the cluster centre has been found between our old and young PMS stars. The younger stars are more concentrated towards the centre of the cluster, while the older stars are typically farther away. This confirms the result from B10, who also reported a difference in the spatial distribution of young and old PMS stars in NGC 3603. Similar bimodal distributions have been found in other studies of PMS stars in massive SFRs (De Marchi et al. 2011b; Gouliermis et al. 2012; De Marchi et al. 2024).

This may simply be a result of older PMS stars having more time to move through the cluster compared to younger PMS stars. The bimodal distribution provides encouraging evidence that there is a meaningful difference between the two samples of PMS stars, and that our age determinations are robust. Our observations are intrinsically biased to not observe stars close to the cluster centre, as the crowding is too intense even for JWST. As such, we may miss older PMS stars located in the cluster centre. Of course, this bias also applies to the younger PMS stars, which may serve to cancel out any relative shift between the two distributions.

7.10 M_acc versus nebular H₂ emission

We have found a relationship between MIacc and nebular H₂ emission. Sources that show enhanced Ṁ_acc relative to M, and age are usually associated with higher than average levels of nebular H2 emission, suggesting an environmental influence on Ṁ_acc. This relationship could be a result of high interstellar gas density that facilitates late-stage infall onto the protoplanetary disk, as has been seen in recent studies (e.g. Cacciapuoti et al. 2024). This infall could replenish the disk, allowing accretion to occur at high levels for longer. Additionally, the higher gas densities may help to protect the disks from external photoevaporation, thus retaining a larger fraction of primordial disk material to accrete.

It is also possible that the relationship works in the opposite direction. Rather than enhanced levels of Ṁ_acc being a result of higher molecular gas densities (inferred from H2 emission), the H2 emission could instead be a result of enhanced levels of turbulence, X-ray and UV emission due to accretion and ejection processes from the stars. However, if the enhanced nebular H2 emission were actually a result of high accretion rates, we would expect to see a similar trend in the nebular H_a emission, but this is not the case.

In W24 the authors found a significant correlation with the mass-normalised accretion rate and gas density traced by 100 am emission. We believe that we have discovered a similar relationship for the young stars in NGC 3603. These early results both indicate that the environment in which a star and protoplane-tary disk form can significantly influence their evolution. While external photoevaporation is being considered more and more as an important environmental factor related to the formation and evolution of protoplanetary disks and planets themselves (e.g. Somigliana et al. 2020; Winter & Haworth 2022; Maucó et al. 2023; Aru et al. 2024), the role of enhanced gas and dust densities associated with massive SFRs has received less attention. External photoevaporation serves to dramatically reduce the lifetime of protoplanetary disks, becoming the dominant disk dispersal mechanism even for modest UV radiation fields (e.g. Miotello et al. 2012). A natural explanation as to how protoplan-etary disks survive for a long time in high UV environments is from protection and/or replenishment from the interstellar medium, which is systematically more dense in massive SFRs (e.g. Heiderman et al. 2010). Regions in which external pho-toevaporation is expected to be important will therefore also be associated with a higher density interstellar medium. These two properties may serve to offset each other, or as our results suggest, the high gas densities appear to dominate, allowing protoplanetary disks to survive for >10 of Myr. The relationship that we have discovered between nebular H₂ emission and M_acc serves as further motivation to understand how these various environmental factors ultimately influence the birth and evolution of stars and planets.

8 Conclusions

We have presented our analysis of 42 PMS stars and their proto-planetary disks in NGC 3603 with JWST NIRSpec spectroscopy. Our main conclusions are the following:

1. The nebular contribution to our stellar spectra has been removed using a scaled subtraction technique that accounts for spatial variations in nebulosity between adjacent micro shutters;

2. 42/100 stellar spectra exhibit strong hydrogen recombination lines after nebular subtraction;

3. All PMS sources have been spectrally classified for the first time. The majority of CTTS have spectral types M, K, and G. The remaining are intermediate-mass stars, including a number of Herbig AeBe type stars;

4. Emission from Fe II and CO bandheads are detected towards a number of young, intermediate mass sources. The presence of these features indicates that, even under the influence of external irradiation, a significant gas-rich inner disk can remain around sources with sufficient mass;

5. No circumstellar H₂ emission is detected towards any source. The lack of this relatively common feature in CTTS and Herbig AeBe spectra may be the result of external irradiation of the outflow that H₂ typically traces;

6. 26/42 sources have Class III SEDs based on their spectral index α. 22 of those source spectra were best fitted with zero veiling emission added to their spectrum. The lack of excess emission indicates that these sources are at late stages of evolution, or their disks have been significantly truncated;

7. The ages, masses, and radii of all 42 sources have been determined with stellar isochrones and evolutionary tracks;

8. The majority of sources have ages ≤5 Myr. 12/42 stars have ages 10-20 Myr. This is significantly longer than the expected disk dispersal timescale of ≤7 Myr;

9. The mass accretion rate of our sources spans five orders of magnitude. There is a strong relationship between M_acc and M_* for stars >1 M_⊙. This relationship is significantly steeper than found in previous studies;

10. M_acc is greater for stars in NGC 3603 for a given age and decreases with age more slowly compared to low-mass SFRs.

11. The old PMS stars (>10 Myrs) appear significantly further away from the cluster centre compared to the young PMS stars;

12. After removing the influence of both age and mass on M_acc, we have found a relationship whereby sources associated with high density molecular gas tend to show enhanced M_acc. This apparent environmental influence on M_acc may help to explain the high accretion rates and old ages of many of our sources.

Data availability

Best fitting model spectra available on Zenodo.

Acknowledgements

We wish to thank the referee for their helpful comments which have improved the quality of this work. Software: SciPy (Virtanen et al. 2020); emcee (Foreman-Mackey et al. 2013); NumPy (Harris et al. 2020); Matplotlib (Hunter & Dale 2007).

Appendix A Spectrum averaging

Each source was observed three times as a result of the nodding pattern employed. We developed an averaging procedure that maximises the S/N of the spectra while also efficiently removing spikes and cosmic rays from the spectra that were missed during data reduction. Our averaging procedure works by comparing and averaging the 2D-rectified and 1D extracted spectra across the three nods. The 2D spectra are 2D arrays with seven pixels in the cross-dispersion direction (rows) and 3817 pixels in the dispersion direction (columns). Our procedure is as follows: The 2D-rectified spectra are averaged together by calculating the median of each pixel in the 2D arrays. The resulting averaged 2D-rectified spectrum is then optimally extracted. We call this spectrum S₁. Next, the procedure goes back to the individual 2D-rectified spectra, and optimally extracts 1D spectra from them. These 1D spectra are then averaged into a single spectrum by finding the median of each wavelength bin. This spectrum we call S₂.

This produces two median spectra: S₁ and S₂. S₁ has typically lower S/N compared to S₂, but virtually all spikes are successfully removed during the averaging process. S₂ on the other hand shows typically higher S/N, but many spikes make it through the averaging process. The reason for S₂ being more resilient to spikes than S₁ is illustrated in Figure A.1. On the right of Figure A.1, which represents the averaging and collapsing process for S₂, we see that two of three 2D spectra contain a spike (red pixel). In this case, the spikes occur in the same detector column for each spectrum and hence at the same wavelength, but in different rows. When these 2D spectra are collapsed, the spikes pass through to the 1D spectrum. When these 1D spectra are averaged, since the spikes occur at the same wavelength, the final averaged spectrum also contains the spike. On the left of Figure A.1, representing the collapsing and averaging process for S₁, since the spikes are present in the same column, but not the same row, they are filtered out during the averaging process and hence are not passed to the collapsed 1D spectrum. For the spikes to pass through in this case, they would need to occur in both the same detector column and row, which is unlikely to occur since these spikes are random, due to cosmic rays and spurious pixels.

Fig. A.1 Left panel, (top) shows three 2D spectra. Two contain a spike due to spurious pixels. (Middle) resulting 2D spectrum after averaging. Spikes are no longer present, having been averaged out. (Bottom) final median 1D spectrum with no spike. Right panel, (top) shows three 2D spectra. Two contain a spike due to a spurious pixel. (Middle) resulting 1D spectra after collapse. Both spikes are passed from the 2D to 1D spectra after collapse. (Bottom) final median 1D spectrum containing the spike.

By comparing S _i and S₂, the spikes in the higher S/N spectrum S₂ can be removed. This procedure was used for both the stellar spectra and the nebular spectra.

Appendix B Fitted Hydrogen Lines

Fig. B.1 Paα emission line profiles of the PMS sources. The best fitting Gaussian profile is shown as a red line. Note that source 523 exhibits a central inversion in its emission line profile, and so no Gaussian fit was performed for this source.

Appendix C Emission line properties

References

All Tables

All Figures

	Fig. 1 Relationship between the EW of He 1 1.869 μm with: left, Paα; centre, Br7; and right, Br6.
In the text

	Fig. 2 Trace plot of source 995.
In the text

	Fig. 3 Corner plot showing the distribution of parameter values.
In the text

	Fig. 4 Best fitting model spectrum shown in red, overplotted on top of the source spectrum shown in black.
In the text

	Fig. 5 Best fitting optical model spectrum for source 152 based on its MUSE spectrum. The black line corresponds to the MUSE source spectrum, the red line to the best fitting model spectrum, while the HST photometry is shown in grey.
In the text

	Fig. 6 As above for source 354.
In the text

	Fig. 7 Histogram of evolutionary classes for our sample.
In the text

	Fig. 8 HRD of the PMS sources using spectroscopically determined Teff and L*. The red points were classified spectroscopically. The orange diamonds were classified based on SED fitting. The lcr uncertainties are represented with red error bars. A selection of isochrones are shown as turquoise dashed lines. The ages are given in Myr on the right hand side.
In the text

	Fig. 9 Location of accreting PMS stars in our sample. Sources with ages <10 Myr in red. Sources with ages >10 Myr in blue.
In the text

	Fig. 10 Results of bootstrapping analysis. The projected distances of young stars in red. Old stars in blue. 95% confidence intervals shown as dashed red and blue lines. The median distance is shown as a dashed black line for each distribution.
In the text

	Fig. 11 Accretion rate as a function of stellar mass. The size of each circle is proportional to age. The largest circle corresponds to 19 Myrs. The best fitting line to our observations is shown in brown, while the best fitting line from Dll is shown in black.
In the text

	Fig. 12 Accretion rate as a function of stellar age. The size of each circle is proportional to M⊙ The largest circle corresponds to 7 MQ. The best fitting line to our observations is shown in brown. The best fitting line from Sicilia-Aguilar et al. (2005) is shown in black.
In the text

	Fig. 13 Normalised accretion rate as a function of the nebular H2 brightness around each source. Sources associated with above average levels of nebular H2 brightness are coloured red.
In the text

	Fig. 14 Relationship between the line flux of Br7 and Paα. The brown line represents the line of best fit and the grey shaded area shows the uncertainty of the fit.
In the text

	Fig. 15 Relationship between Lacc and Paα. The brown line represents the line of best fit and the grey shaded area shows the uncertainty of the tit.
In the text

	Fig. 16 Five PMS sources that display Pfund recombination lines. Their spectra have been normalised to unity, and offset in this figure for visual clarity. The Pfund lines are marked with a green dashed line. Br6 is also present, marked with a blue dashed line. Note that a number of pseudo absorption features are artificially created as a result of our normalisation.
In the text

	Fig. 17 The normalised spectra of the five CO bandhead emission sources. The spectra have been offset for visual clarity.
In the text

Fig. A.1 Left panel, (top) shows three 2D spectra. Two contain a spike due to spurious pixels. (Middle) resulting 2D spectrum after averaging. Spikes are no longer present, having been averaged out. (Bottom) final median 1D spectrum with no spike. Right panel, (top) shows three 2D spectra. Two contain a spike due to a spurious pixel. (Middle) resulting 1D spectra after collapse. Both spikes are passed from the 2D to 1D spectra after collapse. (Bottom) final median 1D spectrum containing the spike.

In the text

	Fig. B.1 Paα emission line profiles of the PMS sources. The best fitting Gaussian profile is shown as a red line. Note that source 523 exhibits a central inversion in its emission line profile, and so no Gaussian fit was performed for this source.
In the text