© The Authors 2025

1. Introduction

Protoplanetary discs around young stars are the site of numerous highly studied astronomical phenomena, including magnetospheric accretion and accelerating disc winds powered by magneto-centrifugal forces, and as the name suggests, they are the formation sites of planets. Understanding the evolution of protoplanetary discs is crucial in both the context of star formation and planet formation. A fundamental push and pull that regulates the lifetime of protoplanetary discs is the balance between mass-loss through winds and mass accretion onto the central star.

Accretion rates of stars were originally measured by fitting shock models to the near-ultraviolet (NUV) spectra of young low-mass stars, known as Classical T Tauri Stars (CTTSs) (Gullbring et al. 1998). The physical picture assumed here is centred around a magnetically driven accretion paradigm in which the star’s strong magnetic field pushes against and truncates the protoplanetary disc out to a few stellar radii (e.g. Hartmann et al. 2016). Material at the disc’s surface is channelled along magnetic field lines from the disc to the central star. During its near free-fall towards the star, the supersonic material shocks the stellar surface and releases its gravitational potential energy as radiation. Initially this is primarily in the form of X-rays, which are quickly absorbed by the stellar surface and re-radiated at longer wavelengths in the NUV and blue-optical. This additional source of NUV light on top of the intrinsic photospheric contribution can be observed and measured, and from this so-called accretion luminosity, the mass accretion rate can be determined. This is the only direct method of measuring the accretion luminosity, and while it unambiguously probes the accretion shock, the reliance on NUV wavelengths makes this approach impractical for sources that are deeply embedded or located at great distances due to high levels of extinction.

An invaluable calibration was discovered by Muzerolle et al. (1998b), who showed that the line luminosity of the relatively bright near-infrared (NIR) hydrogen line Br₇ scales tightly with the accretion luminosity for a sample of 19 CTTSs. With extinction more than ten times lower in the NIR compared to the NUV and without the need for shock modelling, this line luminosity calibration opened the floodgates for the efficient and straightforward estimate of accretion rates for large samples of sources. These calibrations have since been updated and expanded, most notably by Herczeg & Hillenbrand (2008), Alcalá et al. (2014, 2017). We also showed in an earlier work that these calibrations could be extended towards the brightest NIR line – Pa_α (Rogers et al. 2024b), thus making estimating the accretion rate highly accessible in the era of JWST for extremely distant and embedded sources.

All of this work relies on the assumption that accretion is driven magnetically. First theorised by Koenigl (1991) and unambiguously confirmed observationally by Bouvier et al. (2007) for CTTSs, the MA paradigm for Herbig AeBe stars has not reached the same consensus. Unlike CTTSs (spectral types M through G), stars of spectral type A and earlier are not expected to possess a convective envelope, which is thought to be the source of the strong and ordered magnetic fields present around a CTTS (e.g. Kageyama & Sato 1997; Johns-Krull 2007). Despite lingering questions over their origin, magnetic fields have been measured towards dozens of Herbig AeBe stars, albeit with significantly lower field strengths compared to CTTSs (Mendigutía 2020). Other spectroscopic signatures of infalling matter have also been seen towards Herbig AeBe stars, such as red-shifted absorption profiles in He I and H I lines (Cauley & Johns-Krull 2015), though they are less prevalent compared to the occurrence rate in CTTSs. In general, magnetospheres, if present, are expected to be smaller in Herbig Ae stars compared to CTTSs and may be entirely absent for early Be stars (Wichittanakom et al. 2020; Vioque et al. 2022).

Disc winds are ubiquitous around young stars. However, as with accretion, less is known about winds from Herbig AeBe stars compared to CTTSs. Powerful collimated jets and outflows have been spatially resolved for many nearby sources, with large-scale jets being the defining physical feature of the class of objects known as Herbig-Haro objects (Reipurth & Bally 2001). The exact launching mechanism of these winds is still under debate, although there is a growing consensus that these outflows are the result of a magneto-centrifugal wind – also referred to as a magnetohydrodynamical (MHD) wind (e.g. Tabone et al. 2022) – first theorised by Blandford & Payne (1982). These winds are expected to be crucial in the removal of angular momentum from the disc and are theorised to ultimately facilitate accretion by ejecting material, causing outer disc material to migrate inwards and replenish the inner gas disc where accretion takes place (Frank et al. 2014). Studies focusing on disc winds from Herbig AeBe stars, however, have presented contradictory results. In Kraus et al. (2008), the authors used high spatial resolution interferometric observations of the Br₇ line to infer the emission mechanism. They found that four out of five sources were consistent with emission from a stellar or disc wind, with the final source showing more compact Br₇ emission more consistent with emission from the MA flow. In Kreplin et al. (2018), using a similar approach as Kraus et al. (2008), the authors also found that Br₇ emission from the Herbig AeBe star MWC 120 was consistent with a disc wind. In Tambovtseva et al. (2016), the authors used non-local thermodynamic equilibrium modelling of Br₇ emission from a Herbig Ae star and also found that the Br₇ emission likely originates in a disc wind, not an MA flow. In contrast to the above studies, Cauley & Johns-Krull (2014) used the metastable transition of HeI $\lambda 10830\mathrm{\AA }$ as a tracer of mass infall and outflow, which is commonly used for similar purposes with CTTSs (e.g. Dupree et al. 2005; Edwards et al. 2006). In their sample of 56 Herbig AeBe stars, the authors found no evidence of narrow blue-shifted absorption, which traces outflowing material due to disc winds. They did however find red-shifted absorption tracing infalling material for the Herbig Ae stars in their sample, which is indicative of MA. The same red-shifted absorption was not seen towards the Herbig Be stars in their sample. These results do not paint a clear or consistent picture, and they highlight the need for further studies to investigate the accretion and outflow properties of Herbig AeBe stars.

There is no doubt that both winds and accretion give rise to hydrogen emission lines, but there is not yet a consensus on which of the two processes dominates the hydrogen line emission for a given stellar mass and stage of evolution. Using the JWST NIRSpec Micro-Shutter Assembly (MSA) (Ferruit et al. 2022), we obtained 100 stellar spectra, including a number of young intermediate mass stars and Herbig Ae stars. These stars all reside in the giant Galactic star forming region NGC 3603 located 7.2 ± 0.1 kpc away (Drew et al. 2019). We present five of these sources here, which all show rich emission line spectra featuring large unbroken sections of the Brackett and Pfund hydrogen series as well as the strong Pa_α line. We performed a kinematic analysis on the hydrogen lines for each source in order to constrain the physical origin of these lines.

In Sect. 2, we discuss the target selection. In Sect. 3, the data reduction and post-processing steps are explained and the NIRSpec spectra are shown. In Sect. 4, we present the spectral energy distribution (SED) fitting of our sources. In Sect. 5, we describe our methods for fitting and measuring the emission lines. In Sect. 6, the kinematic analysis is presented. In Sect. 7, we interpret our analysis and suggest a physical origin for the emission lines. We also address alternative emission line mechanisms outside magnetic accretion and winds. In Sect. 8, we summarise our findings.

2. Observations and target selection

2.1. NIRSpec observations

The target selection procedure has already been described in Rogers et al. (2024b) and Rogers et al. (2025). Here we summarise it briefly with a specific focus on the five sources presented in this study. 100 stellar spectra were obtained as part of a NIRSpec Guaranteed Time Observations (GTO) programme (PID = 1225, PI: De Marchi, 2022). The aim of this program was to obtain NIRSpec spectra of stars that showed strong line emission due to accretion related activity, as inferred from photometric H_α measurements (see Beccari et al. (2010)). The five targets presented here represent the brightest, most massive and strongest accreting of the sample. These sources had not been studied in detail before, and we did not know prior to these new NIRSpec observations that they constituted a sample of intermediate mass T Tauri stars and Herbig Ae stars (see Sect. 4). As such, we did not set out to obtain spectra of these sources in order to perform a kinematic analysis of their hydrogen lines. Rather, upon inspection of their rich emission line spectra, we were naturally motivated to analyse them further. The SIMBAD names and coordinates of these sources are shown in Table 1.

Our observing strategy consisted of opening three neighbouring micro-shutters for each source in a column of the MSA, forming a ‘mini-slit’. Each source was placed in the central shutter, with the upper and lower micro-shutters observing the surrounding nebula. We employed a nodding pattern, nodding the telescope three times such that the source moved from the central shutter, to the upper shutter and finally the lower shutter. This pattern provided us with three exposures for each source, which could then be averaged together. We found in our analysis that combining exposures to create an average spectrum could increase the measured full width at half maximum (FWHM) of the emission lines by up to 50 km s⁻¹. This may be due to slight differences in the wavelength solution between the three exposures, leading to smearing of the emission lines after averaging. Because of this, we opted to instead work on single exposures. For these five bright sources, this did not present any problems as the marginal boost in S/N due to averaging was not needed. The sources were observed with the grating/filter combination F170LP/G235H, providing us with wavelength coverage from 1.66 − 3.2 μm, with a typical resolving power of ∼4000 (the subject of NIRSpec’s spectral resolution is discussed in Sect. 3.3.

2.2. HST observations

As part of our analysis we also included optical and NIR HST/WFC3 photometry (PID: 11360, PI: O’Connell, 2009). The filters employed from this program were F555W, F625W, F814W, F110W and F160W. A photometric study of the pre-main-sequence population of NGC 3603 was presented in Beccari et al. (2010) using the F555W and F814W filters. Our NIRSpec observations serve as a spectroscopic follow up to those observations. We also used the F435W filter from HST/ACS observations (PID: 10602, PI: Maiz-Apellaniz, 2005). Both programs observed the centre of NGC 3603 with fields of view ranging from 123″ × 139″ for the WFC3 IR channel, 160″ × 160″ for WFC3 UVIS channel, and 202″ × 202″ for the ACS. These two programs were not contemporaneous, taking place almost five years apart. This introduces the potential for photometric variability. Likewise, our NIRSpec spectra were not observed contemporaneously with either of the HST programs. In Rogers et al. (2025), we flux calibrated our NIRSpec spectra using the F160W photometry of each source, after performing a small extrapolation to cover the wavelength gap. The NIRSpec spectra were in excellent agreement with F160W photometry, with typical flux correction factors of 1.02  ±  0.11. This is in agreement with the expected flux accuracy of NIRSpec Gordon et al. (2022), strongly implying that variability does not significantly affect the continuum/broadband photometry for these non-contemporaneous observations. We proceeded under this assumption for the remainder of the analysis.

3. Data reduction

3.1. NIPS

The observations 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 observations which are briefly discussed. 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. The dispersed NIRSpec spectra are curved along the detector. Rectification is performed in order to ‘straighten’ the spectrum. By doing this, the spectrum is resampled onto a uniform wavelength grid (each wavelength bin Δλ is the same for every pixel). The rectified spectrum is a count-rate map consisting of 3817 pixels in the dispersion direction and 7 pixels in the spatial direction. The final data reduction step – extraction, simply collapses the 2D spectrum by summing the 7 pixels along each column.

3.2. Nebular background subtraction

The bright hydrogen emission lines originating from the H II region were subtracted from each of the stellar spectra following the approach outlined in detail in Rogers et al. (2025). This approach makes use of the nebular He I line at (λ = 1.869 μm). The nebular spectrum of each source was scaled and then subtracted, such that the He I line was completely removed from the stellar spectrum. We have estimated that the uncertainty this introduces to the final flux of the hydrogen lines is ∼10%. Figure 1 shows the subtracted NIRSpec spectra of the five sources.

Fig. 1. NIRSpec spectra of our five PMS sources. The spectra have been normalised at 1.66 μm. For better visualisation, they have been offset from each other vertically and cut off after the final Pfund line at ∼2.9 μm. Source names are indicated on the left of the spectra. Hydrogen emission lines are marked with vertical dashed lines; the Paschen series is in orange, the Brackett series is in blue, and the Pfund series is in green. The CO bandheads are also indicated with pink dashed lines. A systematic calibration uncertainty impacts the continuum at ∼1.7 μm

3.3. Spectral resolution of NIRSpec MSA

The often quoted resolving power for NIRSpec in high resolution mode is $\frac{\lambda }{\mathrm{\Delta }\lambda }\sim 2700$. This is based on pre-flight calculations, assuming a fully illuminated micro-shutter, with each resolution element being sampled by 2.2 pixels on the NIRSpec detectors. For point sources however, the resolving power increases by a factor of roughly ∼1.5. To explain this effect, the connection between a micro-shutter and the pixels onto which it projects needs to be understood. The micro-shutters have been designed to project geometrically over two detector pixels, ensuring that for a fully illuminated micro-shutter the line spread function (LSF) is Nyquist sampled (Ferruit et al. 2022). A uniformly illuminated micro-shutter can be imagined as containing an infinite series of point sources. Each of these point sources projects to a slightly different region of the detector. When light is dispersed from each of these point sources, the LSFs of neighbouring points blend together, forming a broader LSF that is composed of many intrinsically narrower LSFs. This blending reduces the fundamental spectral resolution that is achievable with NIRSpec. In the case of a point source, this blending does not occur, as there is only one source of light projecting onto the detector. This means that a higher spectral resolution is achievable when observing point sources compared to extended sources. This has the side effect that, for point sources, the NIRSpec LSF is not Nyquist sampled at any wavelength.

Given that all of our sources are point sources, the actual spectral resolution achieved for our observations needed to be determined if we wished to analyse the kinematics of the emission lines. We have used the JWST-MSAfit software from de Graaff et al. (2024) to do this. This software provides forward modelling and fitting of simulated MSA data. It allows the user to first select a grating and filter combination, and then define a source (extended or point-like) and place it in a specified micro-shutter at a specified position within that micro-shutter. A spectrum with 17 uniformly spaced emission lines is generated from this. The FWHM of the emission lines are measured in order to determine what the resulting resolving power is at each wavelength. As already mentioned, the LSF is not Nyquist sampled by the NIRSpec detector in high resolution mode for point sources. This results in a somewhat jagged FWHM curve (see black line in Figure 5), as some line peaks are centred on one pixel, while others are spread over two. The estimated uncertainty of this approach is expected to be between 10 − 20% (de Graaff et al. 2024). Using this tool we have determined that the actual resolving power for our sources is $\frac{\lambda }{\mathrm{\Delta }\lambda }\sim 4300$, at λ = 1.875 μm. This corresponds to a FWHM of ∼70 ± 14 km s⁻¹, assuming an uncertainty of 20 %.

4. Determining the stellar properties of our sources

In order to draw physical conclusions from the emission lines of our sources, the physical properties of the stars needed to be determined. Given that none of the five sources display photospheric absorption lines in their NIRSpec spectra, traditional spectral classification was not possible. Instead we relied upon SED fitting to determine the underlying photospheric properties and extinction of each source. To do this, we have used archival Hubble Space Telescope (HST) photometric observations in order to perform aperture photometry and measure the optical and infrared fluxes of the sources. The filters employed were F435W (ACS), F555W, F625W, F814W, F110W, and F160W (WFC3). We explored the archives for other photometric observations of these sources, but given the distance and the highly crowded field of NGC 3603, we were forced to rely solely on HST. These observations are the only available with the necessary spatial resolution and sensitivity to observe these sources with high S/N, and without confusion in the beam from other nearby stars. 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). This is especially for true for Herbig stars whose hotter photospheres peak at ∼0.5 μm. We have used a Markov Chain Monte Carlo exploration to fit the photometry of each source, determine the best-fitting stellar parameters, and estimate uncertainties on those parameter values.

4.1. Indirect clues about the spectral type of the continuum sources

There were some immediate clues about the underlying stellar properties of our continuum sources. They are highly luminous, representing the highest S/N sources in the sample, and show relatively flat, or rising SEDs (see Sect. 6.2 of Rogers et al. (2025) for a discussion on the spectral index of the sources). The high luminosities and lack of absorption features in their spectra is naturally explained if these sources are of intermediate or high mass, with relatively hot photospheres. On the other hand, if these sources were actually low-mass CTTS, the high luminosities could possibly be explained as a burst of accretion, with the lack of absorption features in their spectra being a result of exceptionally high levels of excess veiling emission from the protoplanetary disc. Veiling emission is often modelled as a blackbody (e.g. Muzerolle et al. 2003; Cieza et al. 2005; McClure et al. 2013; Antoniucci et al. 2017; Alcalá et al. 2021). We have simulated how much veiling emission would be needed in order to render the strong metal absorption feature Mg I λ1.71 μm completely undetectable for a 4000 K star with [Fe/H] = 0.0 and log(g) = 4.00, at the typical S/N of our continuum sources S/N ∼ 100. We found that to render this line undetectable, a veiling factor r_λ of $\frac{F_{\mathrm{disk}}}{F_{\mathrm{star}}}\ge 20$ was required. Such a high veiling factor also dramatically changes the shape of the SED in the NIR, becoming dominated by the shape of underlying blackbody spectrum. We found that at such high levels of veiling there was no combination of veiling, effective temperature (T_eff), and extinction A(V) that could match our observations. From this, a simpler explanation emerges that these continuum sources simply have hot photospheres that intrinsically lack strong metal absorption lines, and are dominated by hydrogen absorption. These hydrogen absorption lines are entirely filled in by accretion related emission, leaving the spectrum with no detectable absorption lines. We tested this by fitting the SED of our sources in order to determine their T_eff.

4.2. Fitting only optical wavelengths

We initially attempted to fit the optical and NIR observations from HST and JWST in order to constrain the photospheric (T_eff), extinction (A(V)) and disc (T_veil, r_λ) properties simultaneously. We ultimately found that the observations available for these sources were not well suited for such a detailed analysis, and by attempting to fit both stellar and disc features, we introduced too much flexibility for our limited wavelength coverage. We lack observations longwards of ∼3 μm which would constrain the shape of the disc’s SED, and break the degeneracy that we found among our fits. As such, we opted to only fit the optical photometry from HST in order to constrain A(V) and T_eff. By considering only optical wavelengths, we filtered out the majority of emission from the disc, greatly simplifying the analysis (while also returning less physical information). It was still possible to determine r_λ once T_eff was known by comparing the flux of the observed spectrum to the flux of the best-fitting photospheric spectrum. We placed an additional constraint during the fitting procedure that the best-fitting photospheric model must have IR fluxes equal to or below the IR HST photometry fluxes and all NIRSpec fluxes. The flat/rising NIR SEDs of these sources are consistent with moderate to strong excess veiling emission. As such, the model photospheric spectrum should be below the observed spectrum.

We employed a Markov chain Monte Carlo (MCMC) exploration with the python package emcee to determine the best-fitting photospheric model for each source. We have described this procedure for all sources in Rogers et al. (2025). We summarise it here. 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.

We have used the Phoenix stellar models from Husser et al. (2013) in this analysis. We degraded each model to the spectral resolution of NIRSpec, and converted the optical portion of each spectrum into photometric data points by multiplying by the HST filter throughput curves. The 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.

We fixed the metallicity at [Fe/H] = 0.0 and log(g) = 4.00 in order to reduce the number of free parameters and hence the flexibility of the fit. We extinguished our sources using the extinction curve from Gordon et al. (2023) with a total-to-selective extinction ratio of R(V)  =  4.8 (see Rogers et al. (2024a) for a full discussion of the extinction properties of NGC 3603). The parameters that we allowed to vary are T_eff, A(V) and a luminosity scaling term. Given that these sources have not been studied in detail before, we opted for flat (uninformative) priors. We allowed T_eff to vary between 3000 − 12 000 K, A(V) to vary between 0 − 10 mag and the luminosity scaling parameter to vary between 10⁻²³ − 10⁻²².

We achieved good fits for our continuum sources, and as we expected from the indirect clues given above, the best-fitting T_eff values tend to be relatively hot, with a median temperature of 7553 K. The best-fitting parameters for our sources are shown in Table 2. The best-fitting model spectrum is shown for source 238 in Figure 2. The best fits for the remaining sources are shown in Appendix F.

Fig. 2. Top-left panel – SED of source 238 from 0.4–3 μm. NIRSpec spectrum, with uncertainty represented by error bar (black). HST photometry (grey). Phoenix stellar model (red). Bottom-left panel – Normalised NIRSpec spectrum (black). Phoenix stellar model (red). Best-fitting model shows no metal absorption lines, consistent with the NIRSpec spectrum. All hydrogen absorption lines are filled in by emission. We have removed the bump at 1.71 μm for aesthetic reasons. Right panel – Zoom-in of the optical and IR HST photometry.

4.3. Comparing SED fitting to spectroscopic classification

Using SED fitting as a method to determine stellar parameters is not as accurate as a full spectroscopic approach. In order to benchmark the accuracy of the parameters returned from the best-fitting SEDs, we obtained archival optical spectra for two continuum sources, source 152 and source 354 in our sample. These sources are not part of this kinematic analysis, as only a few hydrogen lines are seen in emission. Nonetheless, they have been classified in an identical manner to the five sources being discussed here, and so should indicate how accurate our SED fitting approach is. Their optical spectra were obtained with MUSE by Kuncarayakti et al. (2016). Sources 185 and 251 were also observed during this program, but after inspecting their spectra, we found that all optical hydrogen lines were filled in with emission, just as their NIRSpec spectra are, making spectroscopic classification with these spectra impossible. For sources 152 and 354 however, both H_4 − 2 and a portion of the Paschen series are strongly detected in absorption. We performed fits of their optical spectra with the Phoenix stellar models to constrain T_eff and A(V). Table 3 shows the best-fitting T_eff and A(V) from the SED fitting and spectroscopic fitting. The best-fitting optical spectra for these sources are shown in Appendix C. In both cases we underestimated the true T_eff of the sources, which have optical spectra consistent with spectral type A. This appears consistent with our indirect arguments given above. Exceedingly high levels of veiling are required to remove metal absorption lines from the NIR spectra of cooler stars. So much so that the entire SED shape would become incompatible with our observations. Invoking a high T_eff naturally explains the lack of metal absorption, and is broadly consistent with the relatively high temperatures that were returned by our SED fitting approach. From this, we conclude that sources 185, 238 and 823 are Herbig Ae type stars, while sources 251 and 469 have spectra consistent with spectral type G and F respectively. Given the limitations of SED fitting, it is possible that we have underestimated the temperatures of these sources. It is less likely that any of these five sources are cooler than the best-fitting SED temperatures. Longer wavelength observations from JWST MIRI would greatly help in constraining the disc emission and allow for simultaneous fitting of the photospheric and disc properties. Alternatively, high spectral resolution optical spectroscopy could allow for the wings of filled in hydrogen absorption lines to be detected (such as in Fairlamb et al. 2015), enabling direct spectral classification for these sources.

4.4. Assessment of uncertainties

The uncertainties returned from our MCMC exploration are statistical uncertainties based on 16th–84th posterior distribution percentiles, the Bayesian equivalent to 1σ uncertainties. They represent the range of plausible parameter values given both our observations as well as our priors. Given that we used uninformative priors with a wide range of permitted values for each parameter, the priors do not contribute significantly to the final best-fit values, and ensure that our uncertainties are conservative, rather than unrealistically small. We feel this is reasonable given how little was known about these sources a priori. Our error bars do not reflect systematic uncertainties, such as those inherited from the interpolation of the Phoenix models (see Czekala et al. (2015) for a discussion on uncertainties arising from spectral model interpolation), the intrinsic accuracy of the Phoenix models themselves, or systematic uncertainties related to imperfect calibration of the NIRSpec spectra.

In Rogers et al. (2025) we employed the same MCMC exploration, but fit both the SED and NIRSpec spectrum of sources with clear photospheric absorption features. The typical uncertainties for T_eff and A(V) in those cases was 58 K and 0.2 mag respectively. Comparing this to the error bars returned from the SED fitting of our continuum sources, we find typical uncertainties on T_eff and A(V) of 603 K and 0.33 mag respectively. The error bars assigned to T_eff and A(V) in this study reflect the greater degree of uncertainty associated with SED fitting.

M_* and R_* of each source were determined in Rogers et al. (2025) by fitting evolutionary tracks to the sources on the Hertzsprung Russel diagram (HRD). Likewise, these error bars reflect only the statistical uncertainties on these parameters. M_*, T_eff and L_* of our sample had been estimated before in Beccari et al. (2010) based on a colour-magnitude diagram using the F555W and F814W filters. In that study, a single value of A(V)  =  4.5 mag was applied to all sources to de-extinguish them. The typical uncertainty on M_* was a factor $\sqrt{2}$, corresponding to uncertainties on T_eff of 1000 − 2000 K (private communication, De Marchi.). As we fit each source individually, using more data points and constraints than Beccari et al. (2010), we believe that the best-fitting parameters returned by our MCMC exploration are the most accurate available for these sources currently.

5. Methods – Measuring the recombination lines

In order to measure the hydrogen emission lines in the spectra, we employed Monte Carlo simulations. This allowed for the straightforward propagation of uncertainties that arose from our processing steps and corrections. To do this, we generated 1000 realisations of each emission line, as well as the corresponding absorption lines from the underlying stellar photosphere. The lines were allowed to vary within their uncertainties for each realisation. The absorption lines were varied by generating a new synthetic photospheric spectrum with varying T_eff, A(V) and r_λ based on our results in Sect. 4. In order to reveal the true emission line profiles, the underlying absorption lines needed to be subtracted. This was done by normalising both the observed spectrum and photospheric model spectrum, after correcting for extinction and veiling. Figure 3 shows a section of the normalised spectrum of source 238 after subtraction of the photospheric absorption lines.

Fig. 3. Section of the normalised spectrum of source 238. The grey thin line shows the observed spectrum. The red dashed line shows the photospheric model spectrum after correcting for veiling. The black thick line shows the observed spectrum after subtraction of the photospheric absorption lines.

We performed a fit of a Gaussian profile to each realisation of the emission lines using the non-linear least squares fitting routine curve fit from SciPy. We calculated the equivalent width (EW) and FWHM of the best-fitting Gaussian, and converted the EW to a flux by multiplying it by the adjacent continuum. The final flux and FWHM of each line is the median of the 1000 measurements and we took as its uncertainty the standard deviation of the 1000 measurements. Appendix D lists the emission line EWs and luminosities for each source. In Figure 1, there is a noticeable bump in the spectra at ∼1.71 μm, coincident with the Br₁₀ line. This is a calibration artefact that is introduced after the filter throughput correction step is performed in the pipeline. The line profiles, uncertainties, and best-fitting Gaussian profiles for source 238 are shown in Figure 4. The line profiles for the remaining sources are shown in Appendix E.

Fig. 4. Fitted hydrogen line profiles for source 238. The black solid line is the normalised line profile with the continuum set to zero. The error bars represent the 1σ uncertainties on the flux. The red dashed line is the best-fitting Gaussian profile. The transition name is given in each plot.

6. Results – Line kinematics

6.1. FWHM of hydrogen emission lines

Emission lines are broadened through numerous processes, including natural broadening, thermal broadening, rotational (Doppler) broadening and pressure broadening (eg. Stark broadening). In general, these processes contribute ≤100 km s⁻¹ of broadening (although in some cases it can be significantly more as we discuss in Sect. 7.4). In the circumstellar environment, gas moving at high velocities due to infall and outflow processes can introduce significant line broadening, of order hundreds of km s⁻¹. As such, the FWHM of emission lines from PMS stars is typically interpreted as reflecting the velocity of the gas that produced the emission line. Narrow lines then originate from slow moving gas (e.g. molecular outflows of H₂ and CO with FWHM ≤ 10 km/s⁻¹), while broad lines originate from high-velocity gas (e.g. H I, FWHM ≥ 100 km s⁻¹). Figure 5 shows the FWHM of the Paschen, Brackett and Pfund emission lines from source 238. We have also plotted the resolving power of NIRSpec obtained in Sect. 3.3. Although we focus on source 238 for the remainder of this analysis and discussion, the same kinematic behaviour is seen for all five sources. The FWHM diagrams of the other sources can be seen in the Appendix A.

Fig. 5. Full width at half maximum of each hydrogen emission line of source 238. The Paschen series is shown in orange, and the Brackett series is shown in blue. The dashed line is before photospheric subtraction, and the solid line is after photospheric subtraction. The Pfund series is shown in green. The dashed line is before photospheric subtraction. The solid line is after photospheric subtraction. The expected instrumental resolution is shown in black. The systematic uncertainty of the spectral resolution is in pale grey. Transition names are indicated. For the Pfund transitions, only nup is shown for visual clarity.

A striking relationship is immediately apparent. The FWHM for a given series tends to decrease with increasing wavelength. This is evidently not purely a wavelength effect, as lines from different series with similar wavelengths have drastically different FWHMs. Rather, the trend appears to follow the upper energy level n_up of the transition, with higher order transitions (electrons falling from high n_up) tending to have larger FWHMs. This result suggests that the high excitation lines come from the fastest moving gas, while lower excitation lines come from gas with somewhat lower velocities. The Pfund series displays some significant scatter, as the highest order lines are bunched close together in wavelength and are much weaker than the Brackett or Paschen series. This makes the measurement of these lines challenging in terms of identifying the true level of the continuum, and therefore obtaining reliable fits to the weakest lines. The kinematic trend is dramatic enough that it is still clearly present despite these challenges.

6.2. Line profile optical depths

An additional way to investigate the kinematics of the hydrogen emission lines is to take the ratio of two spectrally resolved line profiles, in order to probe the optical depth of the lines at various velocities. This analysis is akin to that performed by Bunn et al. (1995) and Lumsden et al. (2012). We have used two pairs of lines in this analysis, $\frac{{\mathrm{Br}}_{11}}{{\mathrm{Br}}_6}$ and $\frac{{\mathrm{Pf}}_{17}}{{\mathrm{Pf}}_{11}}$. These lines were chosen so that the optical depth of the Brackett and Pfund series could be compared to each other, as well as to compare lines with significantly different n_up/excitation energies. If the lines are optically thin, the line profile ratio should be consistent with Case B recombination (Baker & Menzel 1938). If the emission is completely optically thick, then the ratio simply becomes a function of wavelength and emitting area:

$$\begin{array}{c}\hfill \frac{I_1}{I_2}={\left(\frac{{\lambda }_2}{{\lambda }_1}\right)}^4\left(\frac{S_1}{S_2}\right),\end{array}$$(1)

where I is the intensity of the emission line, λ is the wavelength of the line, and S is the emitting surface area. This equation is taken from Lumsden et al. (2012), page 1095.

One would expect that the emitting area S should be smaller for the high excitation lines. The degree of excitation in the gas increases closer to the central star. A relatively small portion of gas in the accretion flow or wind is excited enough to produce high excitation lines, while a larger, more extended portion of this gas farther away from the central star is only excited enough to produce lower excitation lines. For simplicity however, we assumed that S₁ = S₂, which provides us with a lower limit to the optically thick ratio for each pair of lines. Line ratios that lie between the completely optically thick and optically thin regimes correspond to only one of the two lines being optically thick. For the optically thin limit for a given ratio, we have used the tables provided in Storey & Hummer (1995). We considered the full range of electron temperatures (500 − 30 000 K) and densities (10² − 10¹² cm⁻³), and adopted the smallest value. This represents the ratio that is closest to optically thick while still being consistent with optically thin emission. In other words, an upper limit. The optically thick limit has been computed based on Eq. (1), and as discussed, represents the lower limit.

Figure 6 shows the optical depth diagrams for $\frac{{\mathrm{Br}}_{11}}{{\mathrm{Br}}_6}$ and $\frac{{\mathrm{Pf}}_{17}}{{\mathrm{Pf}}_{11}}$ for source 238. The optical depth diagrams are shown for the remaining sources in Appendix B. The uncertainty of the line ratios are shown as green error bars, based on the uncertainties of the individual emission line fluxes. Due to the limited number of pixels sampling each line, we interpolated the line profiles so that they could be accurately centred and aligned. As a result, it is important not to over-interpret these diagrams. Small scale variations in the line profile ratio could be artificial, and likely do not reflect real physical changes in the line optical depth over small velocity scales. Only the overall shape of the line profile ratio should be considered, and where it is clearly optically thin, and clearly optically thick. From these optical depth diagrams, we can see that for both the Brackett and Pfund lines, the lowest-velocity components of the lines are consistent with optically thin emission. The wings of the lines, however, quickly become optically thick, with the highest-velocity components of each line showing the most optically thick ratio. Interestingly the Pfund lines show a much larger fraction of the line profile ratio in the optically thick regime compared to the Brackett lines.

Fig. 6. Optical depth diagrams for spectrally resolved emission lines. Top: Line profiles Br11 and Br6 are shown as the dashed and solid black lines, respectively. The line profiles are plotted on a velocity scale, with each line centred at its rest wavelength. Underneath, the green solid line shows the line profile ratio $\frac{{\mathrm{Br}}_{11}}{{\mathrm{Br}}_6}$. The dark red patch represents the fully optically thick values. The medium red patch represents the partially optically thick values. The light red represents the optically thin values. Bottom: Same as above but for Pf17 and Pf11.

7. Discussion – The origin of circumstellar hydrogen emission lines

7.1. Magneto-centrifugal winds

The two most common origins suggested for hydrogen emission lines from young stars are either the accretion flow and/or an accelerating magneto-centrifugal disc wind (Muzerolle et al. 1998a; Kurosawa et al. 2006; Lima et al. 2010). If a magneto-centrifugal wind were the dominant line emission mechanism, the high excitation lines should come preferentially from close to the launching point near the disc surface, where the gas is hottest and most strongly irradiated by the central star. The gas at this point should then have a low velocity, as it has not yet experienced much acceleration. Likewise, one would expect the lower excitation lines to come preferentially from the gas at greater distances from the star, as the irradiation/temperature is still sufficient to produce low excitation lines, and the effective emitting area is significantly larger than at the disc surface. This gas would then have a high velocity, as it has experienced more acceleration. In this picture, the high excitation lines should display narrow profiles, and the low excitation lines should display broad profiles.

This general trend has been observed in Nisini et al. (2004) who observed a very young Class I protostar, known to exhibit powerful jets and outflows, with the lower excitation Br₇ line appearing significantly broader than the higher excitation lines Br₁₂ and Br₁₃. We have analysed archival JWST NIRSpec IFU data of source TMC1A, a class I protostar, with a spatially resolved outflow (see Harsono et al. (2023) for a full discussion of this source). These observations included the same grating and filter combination as our own, allowing for a direct comparison to our observations. Given that TMC1A is spatially resolved, we expect that the resolving power of these observations are intermediate between our own and that of a uniformly illuminated aperture. We have performed the same kinematic analysis on these data, extracting the flux from the datacube using a 3 × 3 pixel aperture. We have made use of only the Paschen and Brackett lines. The Pfund lines in this case are either too weak or located in the strong H₂₀ absorption feature at 3.08 μm, and therefore they could not be measured reliably. Figure 7 shows the resulting FWHMs for each line. The trend seen here is clear. The highest excitation line Br₁₁ has the smallest FWHM, ∼240 km s⁻¹, while the lowest excitation line Pa_α has the highest FWHM, ∼340 km s⁻¹.

Fig. 7. Full width at half maximum of each hydrogen emission line of source TMC1A. The Paschen series is in orange, and the Brackett series is in blue.

In an accelerating wind, the highest-velocity gas far away from the star should have low optical depth, as the density has fallen off by several orders of magnitude compared to the star/disc surface (Lima et al. 2010). Conversely, the gas at the launching point of the wind, with much higher density, should have higher optical depth. The resulting line profile ratios would then show higher optical depth in the line core, becoming optically thinner towards the wings.

In Figure 8, we show the line profile ratio for Br₁₁ and Br₆ for TMC1A. The optical depth diagram of TMC1A is essentially inverted compared to our sources. The line core is the closest to optically thick, while the wings are firmly in the optically thin regime. We note that we have not attempted to correct the spectrum of TCM1A for extinction. Doing so would cause the flux of Br₁₁ to increase by a greater amount than Br₆, pushing the core of the line profile ratio further into the optically thick regime.

Fig. 8. Line profiles of Br11 and Br6, respectively shown as solid and dashed black lines. The green solid line shows the line profile ratio $\frac{{\mathrm{Br}}_{10}}{{\mathrm{Br}}_6}$. Line profile ratio uncertainties are shown as green error bars.

The kinematic trends seen here suggest that the hydrogen emission lines in TMC1A come predominantly from an accelerating wind, which is supported by the spatially resolved outflow seen towards this source, as well as the presence of other wind and outflow tracers such as H₂ and Fe [II] emission. This conclusion was also reached by the authors (Harsono et al. 2023).

7.2. Magnetospheric accretion

The trends seen for TMC1A in terms of kinematics and optical depth are in total contrast to what we observed towards our sources in NGC 3603. Thus, we now consider the accretion scenario as the dominant line emission mechanism for the sources in NGC 3603. In this scenario, the high excitation lines should come preferentially from gas in the accretion flow close to the stellar surface. Here the conditions are preferable for high excitation line emission, as the gas is strongly irradiated by the star and the temperature is significantly higher than at the disc surface (Lima et al. 2010). This gas would move with high velocity, as it has been in near free-fall from the disc surface and hence has experienced significant acceleration. The low excitation lines should come preferentially from gas at the base of the accretion flow on the disc surface. Here the temperature and irradiation are still sufficient to produce low excitation emission lines, and due to the funnel-like geometry of the accretion flow the effective emitting area is larger at the base of the flow compared to just before it hits the stellar surface (e.g. Muzerolle et al. 1998a). This gas would move with a lower velocity as it has yet to experience much acceleration. In this picture, the high excitation lines should display broad profiles, as they come from the highest-velocity gas, while the low excitation lines should display narrow profiles, as they come from the lower velocity gas. This is consistent with the emission line profiles that we have measured across three spectral series, for over 15 hydrogen emission lines for source 238, as well as for the other four sources in this sample.

The fastest moving gas in an accretion flow should also be the most optically thick, as the density of the gas in the accretion flow increases as it approaches the central star. Again, this is consistent with the line profile ratios measured for our sources.

To further test the accretion scenario, we have calculated what the free-fall velocity (v_ff) would be for each source. v_ff is given by the equation

$$\begin{array}{c}\hfill v_{\mathrm{ff}}=\sqrt{2G{M}_{\ast }(\frac{1}{R_{\ast }}-\frac{1}{R_{\mathrm{trunc}}})}cos(i)sin(\theta ),\end{array}$$(2)

where i is the inclination angle and θ is the angle of the accretion material with respect to the equatorial plane. For all calculations we fixed it at $\theta =\frac{\pi }{2}$, corresponding to material reaching the stellar surface at the poles of the star. R_trunc is the truncation radius, the radius at which the magnetosphere reaches into the disc. R_trunc is often assumed to be equal to 5R_* for CTTS (Hartmann et al. 2016), but is likely smaller for Herbig AeBe stars, a result of their weaker magnetic fields (Mendigutía 2020). We have assumed a truncation radius of 3  ±  1R_* to account for this. v_ff is determined by M_* and R_*, which we have determined via SED fitting of our sources. Given the higher degree of uncertainty associated with SED fitting, as well as the systematic uncertainties associated with stellar evolutionary models, we have taken the conservative approach of employing 3σ uncertainties for M_* and R_* in order to determine uncertainties on v_ff.

Table 4 lists v_ff of each source when the gas makes contact with the stellar surface. We have considered three inclination angles of 0° , 30° , and 60°. An observer would measure a maximum line of sight velocity from the accretion flow when viewing the system face-on, with i  =  0°. We also list the FWHM of the Br₁₁ line, which represents a line from some of the highest-velocity gas that we could measure with low statistical uncertainty.

The maximum v_ff are consistent with the FWHM of Br₁₁ for all sources. This demonstrates that gas in a MA flow in free-fall is capable of producing the observed line widths for our sample. This result, in combination with the kinematic trends that we have presented strongly point towards MA over magneto-centrifugal winds. However, other line emission mechanisms can be present around young stars, and so these must be addressed before we can claim that MA is the dominant line emission mechanism for our sources.

7.3. Alternative line emission mechanisms

We have focused on accretion and magneto-centrifugal star/disc winds to try and explain the hydrogen emission lines observed towards our Herbig AeBe stars. There are, however, alternative mechanisms that could potentially explain the kinematic and optical depth results that have been presented so far. They include (1) line driven stellar winds (Kudritzki & Puls 2000), (2) internal/external photoevaporative winds (Williams & Cieza 2011; Winter & Haworth 2022), (3) a hot gaseous inner disc (Muzerolle et al. 2004), and (4) boundary layer accretion (Lynden-Bell & Pringle 1974).

7.3.1. Line-driven stellar wind

Line driven stellar winds are common around massive stars, whose stellar luminosity exceeds the Eddington limit (Abbott 1982). Here stellar photons impart enough momentum to gas on the stellar surface to overcome the gravitational potential of the star, driving a decelerating wind. A decelerating wind could qualitatively match the kinematic results presented. However, although our stars are more massive and hotter than CTTS, their luminosities are still well below the Eddington limit, and hence line driven winds are not expected around our sources.

7.3.2. Internal photoevaporative wind

X-ray and far-ultraviolet radiation from the central star can drive thermal, decelerating winds from the disc surface (Owen et al. 2012). These winds can produce hydrogen emission lines that could in principle display the same kinematic trends that we have measured. Ercolano & Owen (2010) computed theoretical line luminosities for many hydrogen recombination lines forming in a photoevaporative wind. For Br₇, the line luminosity $\frac{L({\mathrm{Br}}_7)}{L_{\odot }}$ ranges from 3.05 × 10⁻⁹ erg s⁻¹ to 1.26 × 10⁻⁷ erg s⁻¹, for X-ray luminosities L_X of 10^29.3 erg s⁻¹ to 10^30.3 erg s⁻¹, respectively. These calculations were carried out for CTTS. X-ray luminosities are generally higher for Herbig AeBe stars compared to CTTS, with L_X up to ∼10³² erg s⁻¹ (Zinnecker & Preibisch 1994; Hamaguchi et al. 2005). This would likely lead to larger internal photoevaporative Br₇ line luminosities. However, these values are still 4 − 6 orders of magnitude lower than the line luminosities measured for Br₇ in our sample. As such, we do not expect that internal photoevaporative winds could dominate the hydrogen line emission for our sources.

7.3.3. External photoevaporative wind

Since our sources are located in a massive star forming region, external photoevaporation could contribute towards the hydrogen line emission. The massive stars of spectral type O and B at the centre of NGC 3603 irradiate our sources, which can result in a thermally driven cocoon of gas and dust escaping the protoplanetary disc, with an ionisation front around the source resulting from extreme-ultraviolet photons. This cocoon morphology is seen commonly around the PMS stars in Orion known as proplyds (e.g. Ricci et al. 2008). At a distance of 7.2 kpc, it is not possible to spatially resolve these sources with current facilities, and so it is not obvious whether they possess this cocoon morphology or ionisation front. However, in the case of Orion, where the detection of these features is unambiguous, the hydrogen emission lines from the ionisation front resemble the nebular emission lines, with narrow line profiles (FWHM ≤ 50 km s⁻¹). Furthermore, our nebular subtraction has removed any contribution from the ionisation front thanks to the He I scaling method, as the He I emission line strength would have a contribution from both the extended nebular gas as well as the ionisation front. By removing the He I emission, we have removed the contribution of both the nebular emission and any possible ionisation front emission. Based on the broad lines that we measure from our sources, and the subtraction method that we have developed and employed, an externally driven photoevaporative wind cannot be the dominant emission line mechanism in this case.

7.3.4. Hot inner gaseous disc

The hot inner gaseous disc may contribute towards the observed hydrogen emission lines. In their seminal work Muzerolle et al. (2004) showed that the inner disc is expected to become an important contributor of emission for mass accretion rates of ≥10⁻⁷ M_⊙ yr⁻¹. In Rogers et al. (2025) we determined Ṁ_acc for all of the sources, including the five presented here, finding accretion rates of 1.6 − 6 × 10⁻⁵ M_⊙ yr⁻¹, well above the level indicated in Muzerolle et al. (2004). The inner disc as an important source of emission is also supported by the presence of CO bandhead emission detected towards three of the five sources (see Figure 1), which likely arises from the warm inner disc (e.g. Ilee et al. 2013). The detection of CO bandhead emission (along with the high values of Ṁ_acc) suggests that there is still ample gas in the inner disc around the stars. In Mendigutía et al. (2017), spectro-interferometric observations of two Herbig AeBe stars revealed spatially unresolved H_α emission, consistent with an origin from a compact inner disc. For one of the two sources, the observations could also be reproduced with a MA model.

If the hydrogen lines were formed in the inner disc, rotational broadening would likely be the dominant line broadening mechanism, since it is the only mechanism capable of producing the FWHMs that we have measured. However, rotational broadening is often associated with a characteristic double-peaked profile in emission lines. Even with the moderate spectral resolution of NIRSpec, double-peaked profiles should be detected, even at low inclinations.

We created a simple toy model of Br₇ emission from the inner disc to demonstrate this. Taking the stellar parameters for source 238 from Sect. 4, we modelled a disc in Keplerian rotation around a star with M_* = 7 M_⊙ and R_* = 12 R_⊙. The disc inner radius R_inner was set to 2R_*, and the outer radius R_outer was set to 1 AU, equivalent to 18R_*, with 100 steps between the two. The Keplerian velocity field was used to calculate the Doppler shift experienced by the line for each step in radius. We modelled the Br₇ emission line as an unresolved Gaussian (FWHM = 50 km s⁻¹). A Doppler shift was applied to the Gaussian for each step in radius and the resulting Gaussians were summed together, producing the final emission line profile, normalised to its peak intensity. The resulting line profile was double peaked. We then downgraded the double peaked profile to the spectral resolution of NIRSpec at the wavelength of Br₇ and also injected Gaussian noise into the line, at the level of 2% (this is the noise level measured in the actual spectrum of source 238). The noisier, lower resolution line profile still exhibits clear double peaked emission line profiles down to inclinations of 20deg. Below this inclination, the double peaked profile disappears and the resulting FWHM reduces to ∼100 km s⁻¹. In Figure 9, the profiles of Br₇ are shown at inclination angles of 10°, 30°, and 70°, respectively. A similar result was also found by Wilson et al. (2022), who computed rotationally broadened emission line profiles for Br₇ for different rates of rotation and inclinations. They showed that for rotational velocities of ∼130 km s⁻¹, which is typical for Herbig AeBe stars (Böhm & Catala 1995), the peak-to-peak separation of double-peaked Br₇ is ∼100 km s⁻¹ for a viewing inclination of 20°, which would be resolved at NIRSpec’s resolution. We note, however, that single peaked emission lines from the Herbig AeBe star HD 100546 have been attributed to emission from a Keplerian disc based on high spectral/spatial resolution observations (Mendigutía et al. 2015), for an inclination angle of i  =  30°. The Br₇ line from that study showed a FWHM of ∼ 200 km/s, similar to the values in our sample. To further test whether the line emission from our sources could find its origin directly from the Keplerian disc, we have calculated the velocity (v_kep) of gas in Keplerian orbit around each source, similarly to our approach in Sect. 7.2. We have compared v_kep with the FWHM of Br₁₁ to test whether the velocities are consistent with the line widths measured for each source. v_kep is defined as

$$\begin{array}{c}\hfill v_{\mathrm{kep}}=\sqrt{G{M}_{\ast }/R_{\ast }}sin(i)cos(\theta ),\end{array}$$(3)

Fig. 9. Line profile of Br7 if it were to originate in a Keplerian rotating inner disc. Inclination angles are 10°, 30°, and 70°.

where G is the gravitational constant, i is the inclination, and θ is the azimuthal angle of the gas in orbit. Here, θ  =  0° corresponds to gas moving towards the observer, and θ  =  180° corresponds to gas moving away the observer. We again considered three values for i of 30° , 60° , and 90°, where i  =  0° corresponds to a face-on disc with no line of sight radial velocity shift and where i  =  90° corresponds to an edge on disc with a maximum line of sight radial velocity shift. We have again employed 3σ uncertainties on M_* and R_* to determine the uncertainties on v_kep. Table 5 shows the maximum v_kep (θ  =  0°) for different inclinations. Sources 185, 238, and 823, the three objects that we have classified as Herbig Ae stars, show FWHMs that are too large to be consistent with emission from the Keplerian disc, even with our conservative approach of employing 3σ uncertainties for the stellar parameters, for a perfectly edge on geometry of i  =  90°, and assuming the disc reaches all the way to the stellar surface. Sources 251 and 469 on the other hand are consistent with emission from both the Keplerian disc and MA.

7.3.5. Boundary layer accretion

An alternative accretion mechanism that has been invoked for Herbig AeBe stars is boundary layer accretion. In this scenario, the magnetic field from the star is either entirely absent or is too weak to truncate the protoplanetary disc. The inner disc connects to the stellar surface via the so-called boundary layer. Disc material must reduce its Keplerian velocity to match the stellar rotational velocity. This reduction in velocity, and hence kinetic energy, is converted into radiation producing the observed accretion luminosity (Mendigutía 2020; Wichittanakom et al. 2020). Given that material must decelerate in order to be accreted, with the highest excitation gas then moving at the slowest velocity, the observed kinematic trends cannot be explained if the bulk of emission comes from the boundary layer. However, in the boundary layer accretion scenario, a Keplerian disc extends close to the stellar surface. It is possible that the bulk of the line emission may come from the disc itself due to its larger emitting surface compared to the boundary layer. As such, even if boundary layer accretion were taking place for our sources, the kinematic information retrieved from the hydrogen emission lines may still be largely governed by the Keplerian disc. Because of this, we cannot rule out boundary layer accretion as a possibility for sources 251 and 469.

7.4. Stark broadening

Up to now, we have assumed that the dominant line broadening mechanism is Doppler broadening, due to the high-velocity bulk motion of the gas. This is often invoked to explain the line widths for CTTS and Herbig AeBe stars, as the velocities from infalling and outflowing material match the line widths well (e.g. Muzerolle et al. 2004). In some cases, line widths have been measured that are too broad to be consistent with Doppler broadening. One such example comes from Muzerolle et al. (1998a), who suggest that Stark broadening could explain the high-velocity wings observed towards H_α, extending out as far as ±500 km s⁻¹. A feature of Stark broadening that is of particular importance in the context of our study is that transitions with high n_up are more susceptible to Stark broadening, due to the weaker Coulomb force between those electrons with the nucleus. As a result, Stark broadening can induce qualitatively similar FWHM trends as those that we have presented. An example of this can be found in Feldman & Doschek (1977), who show the Stark effect on high order Balmer lines from the solar corona. Stark broadening is an important line broadening mechanism in Herbig AeBe photospheres due to their high temperatures and relatively high densities, both of which serve to amplify the Stark effect (e.g. Popović et al. 2001; Dimitrijevic 2007). In Sect. 5, we subtracted the Stark broadened absorption lines from the circumstellar emission lines. We also measured the FWHM of the emission lines before subtraction, in order to understand the impact of the underlying absorption lines. In general, the subtraction led to an increase in the FWHM of the emission lines, but did not change the overall trend seen in Figure 5. This assured that the FWHM trends of the emission lines were not strongly influenced by the absorption component. However, given the potential degeneracy between Stark broadening and doppler broadening in the MA flow, it was prudent to understand whether Stark broadening may be an important broadening mechanism in the circumstellar environment. In Bunn et al. (1995) and Lumsden et al. (2012), the authors argue that Stark broadening may be important for Brackett lines in the circumstellar environment of intermediate and high mass stars, due to the presence of high-velocity wings reaching out to ∼500 km s⁻¹. None of the emission lines in our sample exhibit such broad wings, with the broadest lines having wings out to ∼300 km s⁻¹. In Muzerolle et al. (1998a), the authors found that Stark broadening can impact H_α significantly for high densities n_H ≥ 10¹² cm⁻³ and temperatures T ∼ 10 000 K, but the effect dropped off quickly for higher order lines, with H_β experiencing 50% less broadening compared to H_α. The same conclusion was found by Tambovtseva et al. (2014), who modelled the Br₇ emission line as well as H_α. They reasoned that the high optical depth of H_α leads to significant Stark broadening, while the lower optical depths of Brackett emission lines results in negligible Stark broadening compared to doppler broadening. Finally, Wilson et al. (2022) computed Stark broadened emission line profiles for H_α, Pa₅, Pa₆ and Br₇ and found again that while the effect was significant for H_α, the NIR lines were negligibly affected. Based on these previous results, as well as the lack of very high-velocity line wings in our observations, we do not expect that Stark broadening contributes significantly to our line widths.

We conclude that the line emission from sources 185, 238 and 823 is most consistent with MA. The emission lines for sources 251 and 469 are consistent with both MA and inner disc emission.

7.5. The MA paradigm in Herbig AeBe stars

The observed Paschen, Brackett, and Pfund lines originating from the MA flow for sources 185, 238 and 823 is a somewhat unexpected result. It is surprising that disc wind signatures are not consistent with our kinematic analysis, despite their expected importance for Herbig AeBe stars.

Tambovtseva et al. (2016) performed non-local thermodynamic equilibrium modelling to investigate whether accretion or winds were the dominant emission mechanism for Br₇ from Herbig AeBe stars, concluding that disc winds do indeed dominate. This was supported by a subsequent spectro-interferometric observational study of MWC 120 by Kreplin et al. (2018), who concluded again that Br₇ is formed predominantly in an extended wind around the star, not the compact magnetosphere. Kraus et al. (2008) employed spectro-interferometry to try and spatially resolve the emission region of Br₇. They found that of the five sources observed, only one, source HD98922, displayed compact Br₇ emission, consistent with an origin from a MA flow. It is therefore rather surprising and puzzling that our kinematic results are consistent with MA. One perspective that we have not considered yet is that our sources are all found in a massive star forming region, residing near a population of more than 70 O stars (Moffat et al. 1994; Drissen et al. 1995). This type of environment is significantly different to where the most well studied, nearby Herbig AeBe stars reside. The harsh UV interstellar radiation field may influence the stars’ ability to drive winds/outflows.

In Reipurth et al. (1998), the authors observed four Herbig-Haro (HH) objects close to the Horsehead nebula in Orion. All four of these stars lie in the vicinity of the O9.5 star σ Orionis. In this study the authors found that the typically molecular outflows associated with HH objects are in fact ionised, which they reasoned is a result of external ionisation from σ Orionis. They suggested that jet launching 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 disc itself. Additionally, studies towards the Carina nebula also found evidence of externally irradiated jets (Reiter & Smith 2013; Reiter et al. 2016, 2022) due to the presence of nearby OB stars. Very recent spectroscopic follow up observations have revealed that outflows exposed to external UV irradiation become dissociated, ionised, and eventually totally ablated (Reiter et al. 2024). The massive star content and UV radiation field of NGC 3603 is several orders of magnitude more extreme compared to Orion, and comparable to the Carina nebula. We speculate that the incompatibility between the hydrogen line kinematics and a magneto-centrifugal wind may be a result of the extreme environment that our sources find themselves in. The JWST NIRCam imaging could help answer this question. Outflows with projected sizes of > 400 AU would be spatially resolved, allowing for an assessment of the impact of the external UV field on the formation, morphology and lifetime of jets and outflows in NGC 3603.

The origin of hydrogen emission lines from Herbig AeBe stars is not definitively known, and likely depends on multiple parameters including the age, mass, the accretion rate of the central source, and possibly even their external environment. We are not aware of any other studies in the literature that have investigated the kinematic properties of hydrogen emission lines systematically as we have presented here, and suggest that this approach could be helpful in distinguishing between different emission line mechanisms.

We also wish to emphasise that this kinematic analysis does not rely on having JWST NIRSpec observations. Although NIRSpec is an ideal facility for this type of study given its unbroken wavelength coverage providing many hydrogen lines with high S/N and sufficient spectral resolution, there are many existing archival Herbig AeBe spectra from instruments such as x-shooter that are well suited for this type of study. Some potential challenges with spectra from ground based facilities include gaps in wavelength coverage due to telluric absorption in the NIR, larger extinction correction uncertainties at optical wavelengths, and the importance of other line broadening mechanisms for certain Balmer lines making the kinematic analysis less straightforward. Despite these obstacles, the large sample of archival spectra provides an opportunity to expand this study with a statistically significant sample.

8. Conclusions

We have presented JWST NIRSpec spectra of five intermediate mass sources located in the massive star forming region NGC 3603. The spectra of our sources exhibit many recombination lines (mostly from H I) and lack any detectable absorption lines from their underlying photospheres. We fit their optical photometry with an MCMC exploration to estimate their stellar properties. We have classified sources 185, 238, and 823 as Herbig Ae stars based on the best-fitting T_eff. Sources 251 and 469 have spectra consistent with spectral types G and F, respectively. Sources 185, 238, and 251 exhibit CO bandhead emission. We performed a kinematic analysis on the multiple series of hydrogen emission lines to try to constrain where the line emission originates. Based on the FWHM and optical depth of the spectrally resolved emission lines and considering the typical velocities expected for a number line emission scenarios, we favour an origin in the MA flow for sources 185, 238, and 823. We cannot rule out the possibility of emission from the Keplerian disc for sources 251 and 469. Surprisingly, none of our sources are consistent with emission from a magneto-centrifugal wind, despite this being a commonly invoked mechanism to explain hydrogen line emission for Herbig AeBe stars. This result provides new observational support to the existence of a magnetically driven accretion mechanism around Herbig Ae stars, which can dominate the line emission from these sources. With a small sample size of five, it is unwise to draw generalised conclusions from these sources alone. The kinematic analysis that we have presented and described here can be extended to the large sample of existing archival Herbig AeBe spectra, provided enough hydrogen lines are present with sufficient S/N and spectral resolution. Increasing the sample size would provide insight into the extent to which MA dominates the line emission for Herbig AeBe stars.

Acknowledgments

We wish to thank the referee for their careful review of this work. Their comments and suggestions significantly improved the quality of our analysis. Software: SciPy (Virtanen et al. 2020); NumPy (Harris et al. 2020); Matplotlib (Hunter & Dale 2007); MSAfit (de Graaff et al. 2024)

References

Appendix A: FWHM diagrams

Fig. A.1. Full width at half maximum of each hydrogen emission line. Top left - Source 185. Bottom left - Source 251. Top right - Source 469. Bottom right - Source 823.

Appendix B: Optical depth diagrams

Fig. B.1. Optical depth diagrams of sources 185, 251, 469, and 823. In the case of source 469, Pf13 was used, as neither Pf12 or Pf11 were measured. In the case of source 823, Pf12 was used, as Pf11 was not measured.

Appendix C: MUSE spectra best fits

Fig. C.1. Best-fitting optical model spectrum for source 152 based on its MUSE spectrum. In black is the MUSE source spectrum. In red is the best-fitting model spectrum. In grey is the HST photometry.

Fig. C.2. As before but for source 354.

Appendix D: Emission line properties

Appendix E: Line profiles

Fig. E.1. Line profiles for source 185.

Fig. E.2. Line profiles for source 251. We note that Pf22 − 5 and Pf21 − 5 have been excluded due to low S/N.

Fig. E.3. Line profiles for source 469. We note Pf22, Pf21, and Pf20 have been excluded due to low S/N. Due to the location of the spectrum of source 469 on the NIRSpec detectors, Br8, Pf12, and Pf11 were not measured on the detector.

Fig. E.4. Line profiles for source 823. We note Pf22, Pf21, and Pf20 have been excluded due to low S/N, and Pf11 was not measured on the detector.

Appendix F: MCMC best-fit plots

Fig. F.1. Best-fitting spectrum for source 185.

Fig. F.2. Best-fitting spectrum for source 251.

Fig. F.3. Best-fitting spectrum for source 469.

Fig. F.4. Best-fitting spectrum for source 823.

All Tables

All Figures

Fig. 1. NIRSpec spectra of our five PMS sources. The spectra have been normalised at 1.66 μm. For better visualisation, they have been offset from each other vertically and cut off after the final Pfund line at ∼2.9 μm. Source names are indicated on the left of the spectra. Hydrogen emission lines are marked with vertical dashed lines; the Paschen series is in orange, the Brackett series is in blue, and the Pfund series is in green. The CO bandheads are also indicated with pink dashed lines. A systematic calibration uncertainty impacts the continuum at ∼1.7 μm

In the text

Fig. 2. Top-left panel – SED of source 238 from 0.4–3 μm. NIRSpec spectrum, with uncertainty represented by error bar (black). HST photometry (grey). Phoenix stellar model (red). Bottom-left panel – Normalised NIRSpec spectrum (black). Phoenix stellar model (red). Best-fitting model shows no metal absorption lines, consistent with the NIRSpec spectrum. All hydrogen absorption lines are filled in by emission. We have removed the bump at 1.71 μm for aesthetic reasons. Right panel – Zoom-in of the optical and IR HST photometry.

In the text

	Fig. 3. Section of the normalised spectrum of source 238. The grey thin line shows the observed spectrum. The red dashed line shows the photospheric model spectrum after correcting for veiling. The black thick line shows the observed spectrum after subtraction of the photospheric absorption lines.
In the text

	Fig. 4. Fitted hydrogen line profiles for source 238. The black solid line is the normalised line profile with the continuum set to zero. The error bars represent the 1σ uncertainties on the flux. The red dashed line is the best-fitting Gaussian profile. The transition name is given in each plot.
In the text

Fig. 5. Full width at half maximum of each hydrogen emission line of source 238. The Paschen series is shown in orange, and the Brackett series is shown in blue. The dashed line is before photospheric subtraction, and the solid line is after photospheric subtraction. The Pfund series is shown in green. The dashed line is before photospheric subtraction. The solid line is after photospheric subtraction. The expected instrumental resolution is shown in black. The systematic uncertainty of the spectral resolution is in pale grey. Transition names are indicated. For the Pfund transitions, only nup is shown for visual clarity.

In the text

Fig. 6. Optical depth diagrams for spectrally resolved emission lines. Top: Line profiles Br11 and Br6 are shown as the dashed and solid black lines, respectively. The line profiles are plotted on a velocity scale, with each line centred at its rest wavelength. Underneath, the green solid line shows the line profile ratio $\frac{{\mathrm{Br}}_{11}}{{\mathrm{Br}}_6}$. The dark red patch represents the fully optically thick values. The medium red patch represents the partially optically thick values. The light red represents the optically thin values. Bottom: Same as above but for Pf17 and Pf11.

In the text

	Fig. 7. Full width at half maximum of each hydrogen emission line of source TMC1A. The Paschen series is in orange, and the Brackett series is in blue.
In the text

	Fig. 8. Line profiles of Br11 and Br6, respectively shown as solid and dashed black lines. The green solid line shows the line profile ratio $\frac{{\mathrm{Br}}_{10}}{{\mathrm{Br}}_6}$. Line profile ratio uncertainties are shown as green error bars.
In the text

	Fig. 9. Line profile of Br7 if it were to originate in a Keplerian rotating inner disc. Inclination angles are 10°, 30°, and 70°.
In the text

	Fig. A.1. Full width at half maximum of each hydrogen emission line. Top left - Source 185. Bottom left - Source 251. Top right - Source 469. Bottom right - Source 823.
In the text

	Fig. B.1. Optical depth diagrams of sources 185, 251, 469, and 823. In the case of source 469, Pf13 was used, as neither Pf12 or Pf11 were measured. In the case of source 823, Pf12 was used, as Pf11 was not measured.
In the text

	Fig. C.1. Best-fitting optical model spectrum for source 152 based on its MUSE spectrum. In black is the MUSE source spectrum. In red is the best-fitting model spectrum. In grey is the HST photometry.
In the text

	Fig. C.2. As before but for source 354.
In the text

	Fig. E.1. Line profiles for source 185.
In the text

	Fig. E.2. Line profiles for source 251. We note that Pf22 − 5 and Pf21 − 5 have been excluded due to low S/N.
In the text

	Fig. E.3. Line profiles for source 469. We note Pf22, Pf21, and Pf20 have been excluded due to low S/N. Due to the location of the spectrum of source 469 on the NIRSpec detectors, Br8, Pf12, and Pf11 were not measured on the detector.
In the text

	Fig. E.4. Line profiles for source 823. We note Pf22, Pf21, and Pf20 have been excluded due to low S/N, and Pf11 was not measured on the detector.
In the text

	Fig. F.1. Best-fitting spectrum for source 185.
In the text

	Fig. F.2. Best-fitting spectrum for source 251.
In the text

	Fig. F.3. Best-fitting spectrum for source 469.
In the text

	Fig. F.4. Best-fitting spectrum for source 823.
In the text