© The Authors 2024

1 Introduction

Determining the age, distance, and star formation history of a stellar cluster relies upon an accurate knowledge of the extinction towards the region. Extinction is often characterised with the term $R(V)\equiv \frac{A_{\text{v}}}{E(B-V)}$, known as the total-to-selective extinction. Most regions in the Milky Way follow a typical Galactic extinction law with a total-to-selective extinction ratio R(V) = 3.1. However, young regions with active ongoing star formation tend to exhibit higher values of R(V), implying that the typical size of dust grains is large in these regions compared to the Galactic average. Indeed, many studies towards Galactic HII regions, or starburst galaxies, find systematically higher values of R(V) (e.g. Neckel & Chini 1981; Bautista et al. 1995; De Marchi & Panagia 2014; Pandey et al. 2000; Calzetti et al. 2000) compared to the Galactic average of 3.1. This could be a result of photoevaporation of smaller grains due to radiation from massive stars and/or the coagulation and growth from small to large grains, aided by the higher densities in star forming regions (Swamy 1965). Gall et al. (2014) show that supernovae explosions are also highly efficient means of producing large dust grains.

NGC 3603 is a giant Galactic HII region, 1–3 Myr old, located 7 ± 1 kpc away (Pandey et al. 2000; Melena et al. 2008; Pang et al. 2016). The region hosts an extremely compact OB cluster at its centre, which has made studying the central region in detail difficult due to limited angular resolution. With about one-third the mass in stars of 30 Doradus in the Large Magellanic Cloud (LMC), NGC 3603 is the closest system in the Galaxy to a genuine starburst region Stolte et al. (2006). It is also the optically brightest HII region in the galaxy.

Most studies on the star formation history of NGC 3603 rely upon optical photometry, and make use of either the colour-difference method or two colour diagrams (TCDs) to determine the extinction (e.g. Moffat 1983; Melnick et al. 1989; Pandey et al. 2000; Sung & Bessell 2004), referred to as Mof83, Mel89, Pan00, and Sun04 going forward.

An alternative approach to determine the extinction law and colour excess towards a region is to study the recombination lines of nebular spectra. By comparing the ratios of hydrogen emission lines present in nebular spectra to their theoretical values based on Case B recombination theory Osterbrock & Ferland (2006), the extinction law and colour excess can be determined, based on how much the measured ratios deviate from the theoretical ratios. In this study, we obtained 600 nebular spectra using the multi-object spectroscopy (MOS) mode of NIRSpec on board the James Webb Space Telescope (JWST), using the G235H grating and F170LP filter. This corresponds to a resolving power of $\frac{\lambda }{\text{Δ}\lambda }=2700$, and a wavelength range of 1.6–3.2 µm. MOS mode utilises the novel Micro-Shutter Assembly (MSA), which consists of ~250 000 individual micro-shutters that can be configured open or closed depending on the positions of sources in the field of view.

In Sect. 2 we discuss the target selection. In Sect. 3, the data reduction process is explained. In Sect. 4, the results of our extinction analysis are presented. We also present a high signal to noise ratio (S/N) stacked nebular spectrum in Sect. 4.6. In Sect. 5, we discuss our results, assumptions, uncertainties, and biases, and compare our results to the literature. In Sect. 6 our findings are summarised.

Fig. 1 Region of NGC 3603 with the projected micro-shutter’s angle and dimensions. This RGB image is composed of High Acuity Wide-field K-band Imager (HAWK-I) F212n H2 (Red), HST F160W (Green), and HST F168n (Blue). The stellar sources not discussed in this paper lie within the central micro-shutters, while the upper and lower shutters observe the nebula. The zoomed in section containing the source 3796 illustrates the nodding pattern utilised. JWST nods up and down, moving three neighbouring open micro-shutters north-east and south-west, observing 6 nebular regions in total. The position of the OB cluster is indicated by the blue star and arrows.

2 Selection of targets and observation strategy

The principle aim of the NGC 3603 NIRSpec observations was to gain insight into the formation processes of the pre-main sequence (PMS) stars in the region. One hundred stellar sources were selected for observation with NIRSpec from an initial list of about 10000, originally observed with Hubble Space Telescope’s (HST) Wide Field Camera 3 (WFC3), see Beccari et al. (2010). The 100 sources observed with JWST had been classified photometrically as either PMS (60/100) or main sequence (MS; 40/100), based on Hα excess emission. For our JWST observations, each source was placed within a micro-shutter, with the shutters directly above and below also open, forming a mini-slit of three micro-shutters. The upper and lower shutters were opened in order to measure the nebular background, which is crucial in a region like NGC 3603 where the nebular emission is bright. A three nod strategy was employed, moving the source from the central shutter, to the upper shutter, and finally to the lower shutter so that source spectra could be measured on different regions of the detector.

A major constraint on the final source selection was the placement of each source in a micro-shutter. Sources needed to be relatively isolated, so that contaminant stars were not observed in the neighbouring background shutters. Here a contaminant star is simply any star that happens to fall within one of the open micro-shutters, but is not a target of observation. This phenomenon was impossible to completely avoid in a region as dense with stars as NGC 3603. The impact of contaminants on our nebular spectra is discussed Appendix A.

Figure 1 shows a section of NGC 3603. The bright southeast pillar is visible, along with many projected micro-shutters. Note the sharp drop in nebulosity in the upper left hand corner of the image. The final selection of 100 sources from the original catalogue of 10 000 was primarily motivated by how isolated the sources were, so that contamination could be kept to a minimum. The final observation strategy consisted of three MSA configurations, with three nods per configuration. A full discussion of the MS and PMS stars will follow in Rogers et al. (in prep.).

For this study, we focus on the nebular spectra that were obtained. Given the three nod strategy and two micro-shutters observing the nebular background per stellar source, six nebular spectra were obtained per source, totalling 600 nebular spectra in total.

Due to the nodding pattern of the telescope, certain nebular regions were observed multiple times, producing repeat spectra. The pointing accuracy of JWST is 0.1″ (Lallo & Hartig 2022). Given the distance to NGC 3603 of 7 kpc, this small pointing offset results in sampling regions that are close in angular distance but relatively far apart (~700 AU) in terms of linear distance. For this reason we treat the repeat observations as independent spectra. The typical S/N in the continuum of each spectrum was ~10. The recombination lines are quite strong, having a typical S/N of ≥30, with the strongest lines ≥100.

Fig. 2 Example of four nebular spectra exhibiting a raised continuum level around Br10−4 as a result of the NIRSpec calibration. This effect produces unphysical recombination ratios with respect to Case B recombination, and so this line has not been included in the remainder of the analysis.

3 Data reduction–NIPS

The uncalibrated data was largely reduced with the ESA Instrument Team’s pipeline known as NIRSpec Instrument Pipeline Software (NIPS; Dorner et al. 2016). 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 and spectral extraction, with the final product being the 1D extracted spectrum. One of the final steps before extraction is the rectification of the spectrum. Since NIRSpec spectra do not run perfectly parallel with the dispersion direction, as well as bending slightly along the detector, a rectification step is performed on the so-called irregular 2D spectrum. This step resamples the 2D spectrum on a linear and uniform grid. The rectified spectrum, also called the regular 2D spectrum, is a count-rate map consisting of 3817 pixels in the dispersion direction and seven pixels in the cross-dispersion direction. The extraction step then simply collapses the 2D spectrum by summing the seven pixels along each detector column. During the data reduction process, it became apparent that contamination of the nebular spectra had occurred. This rendered the affected nebular spectra unusable. Contamination reduced the number of usable nebular spectra from 600 to 471. A full discussion of this effect is given in Appendix A.

4 Extinction and colour excess

4.1 Measured and theoretical line ratios

It is possible to measure the extinction in a region via the recombination lines found in the nebular background spectra, (e.g. Calzetti 1997b; Calzetti et al. 1996). Interstellar extinction reddens spectra, extinguishing blue wavelengths more than red wavelengths. The expected intensity ratio of two recombination lines for a given electron temperature T_e and electron density N_e is known from Case B recombination theory. The “B” signifies that the environment is optically thick to all Lyman series photons, that is, these photons are absorbed. The alternative scenario is Case A, where these photons are not absorbed Osterbrock & Ferland (2006). Case B is expected to be valid for HII regions where the density of hydrogen is relatively high (≥100 cm⁻³). The deviation that the measured recombination line ratios exhibit compared to Case B is a measure of the extinction that is taking place.

4.2 The choice of recombination lines

A strong Brackett series line in our nebular spectra could not be used in this extinction analysis − Br₁₀. This was because it is located in a region of the spectrum that is visibly affected by the calibration, and hence we could not measure a reliable line flux. Figure 2 illustrates this, showing four nebular spectra from neighbouring micro-shutters.

The continuum below the line is visibly raised, as indicated with the dashed black line tracing the correct continuum. The absolute change in the flux of Br₁₀ is a few %. However, even a small change in the flux of the line can result in a significant change in the line ratio, which can significantly affect the resulting value of E(B – V). As such, we have excluded Br₁₀ from this analysis.

All of the lines used in this analysis belong to the Brackett series. Using lines from the same series helps to ensure that the nebular conditions (T_e and N_e) being probed are similar. Ideally all of the lines that we have used would have the same upper energy level, so that the excitation potential was the same, meaning that they arise from identical environments. This was not possible from the lines available to us, and so we opted to rely on a single hydrogen series, with all lines relatively close together in wavelength, keeping the physical conditions producing these lines similar. A recombination line list is given in Table 1.

Fig. 3 Five Brackett series line profiles are shown from a representative spectrum, plotted in log-scale. The best fitting Gaussian profile is shown as a red dashed line.

4.2.1 Using Pfund recombination lines to probe extinction

There is a rich forest of Pfund recombination lines present in the nebular spectra. For a number of reasons, it has not been possible to study these lines to determine the extinction characteristics of NGC 3603. These lines lie at longer wavelengths than all of the Brackett lines except for Br₆. This means the extinction experienced by these lines is weaker, and hence harder to measure. The lines are also relatively close together in wavelength, making the differential extinction harder to measure. The lines are also weaker, resulting in larger uncertainties, which also hinders our determination of E(B – V). As we shown in Sect. 4.5.1, due to the overall weak differential extinction that we determine, the Pfund lines are not well suited to measuring E(B – V) in NGC 3603, though we acknowledge that it would provide a very interesting secondary probe of extinction towards the region.

4.2.2 Measuring the recombination lines

For each extracted nebular spectrum, we used Monte Carlo to estimate the line flux and uncertainty for the Brackett recombination lines, similar to that of Ashraf et al. (2023); Riedel et al. (2017). This approach involved fitting the recombination lines many times, allowing the lines to vary within their uncertainties for each iteration. For this study we found that the recombination lines uncertainties were well sampled after 1000 iterations. The flux of each line was determined by calculating the median of the 1000 flux measurements, and the uncertainty of the flux was the standard deviation of the 1000 measurements. We fit the recombination lines with Gaussian profiles using the non-linear least squares fitting routine “curve-fit” from SciPy. In combination with the analytical expression for a Gaussian, this fitting routine returns the best fitting values of the amplitude, standard deviation, and central wavelength of the Gaussian. We determined the equivalent width (EW) of the best fitting Gaussian profile, and converted the EW to a flux by multiplying by the level of the adjacent continuum. In some cases the hydrogen recombination lines were partially blended with another, weaker line. To estimate the impact of these weaker blended lines, we fit a double Gaussian to both lines and found that the relative flux contribution from the weaker lines was typically <1%, so we do not consider it to significantly affect our results. An example of the line profiles of the Brackett series is shown in Fig. 3.

The flux of each hydrogen line was divided by the flux of Br₆, yielding a list of line ratios. The uncertainty of each ratio was found by propagating the error of each flux using the standard error propagation formula. The Br6 line was chosen as the anchor (the line that we compare the other lines to) because it is present in all of the spectra (in some cases, Br₇ is located in the NIRSpec detector gap), and is always strongly detected, resulting in low statistical uncertainties. The measured line ratios then needed to be compared to the expected ratios from Case B recombination. To do this, the line ratios of Storey & Hummer (1995) were used for values of T_e = 10 000−12 500 K and n_e = 100–10 000 cm⁻³. The line ratios are not particularly sensitive to T_e or N_e and so our specific choice of these parameters is not extremely important. Our choice of T_e comes from Esteban et al. (2004); McLeod et al. (2016), who both find typical electron temperatures of about 11 000 K in NGC 3603.

Equation (1) gives the relationship between the colour excess, the line ratios, and the extinction curve, $$E(B-V)=\frac{-\log \left(r_{\text{obs }}/r_{\text{int }}\right)}{0.4\left[\kappa \left({\lambda }_{\text{line1 }}\right)-\kappa \left({\lambda }_{\text{line }2}\right)\right]}\text{mag},$$(1)

where r_obs is the measured ratio of two hydrogen recombination lines and r_int is the theoretical ratio of two hydrogen recombination lines based on Case B recombination. We refer to the term: log(r_obs/r_int) as the “decrement” moving forward. The notation of formula 1 is taken directly from Calzetti (1997a), and may lead to some confusion as the letter r has been used to represent the line ratios, while we use it to represent the extinction curve in this study. In order to minimise confusion, we use lower case to differentiate between r_obs/r_int – the line ratios, and R(V)–the extinction curve. The terms κ(λ_line1) and κ(λ_line2) are equivalent to R(λ_line1) and R(λ_line2), that is, the values of the assumed extinction curve at the wavelengths of the recombination lines in question. We refer to this term as the “extinction curve differences” moving forward. Having numerous recombination lines meant that each set of ratios could provide an independent estimate of E(B – V). The uncertainties of the line fluxes derived from the Monte Carlo iteration method as well as the uncertainties in T_e and N_e were propagated through to each value of E(B – V) using the standard error propagation formulae.

4.3 The choice of extinction curve

If an inappropriate extinction curve is used in Eq. (1) to determine E(B – V), different values of E(B – V) will result from different line ratios. The appropriate extinction curve on the other hand, should produce the same value of E(B – V) regardless of which line ratio is used, though this may not hold for lines with very large wavelength separations. This caveat is fully discussed in Sect. 5. The form of the NIR extinction curve that we have adopted is that of Fitzpatrick & Massa (2009, hereafter FM09), $$\frac{E(\lambda -V)}{E(B-V)}=[0.349+2.087R(V)]\times \frac{1}{1+{(\lambda /0.507)}^{\alpha }}-R(V),$$(2)

where α is the exponent that defines the curvature of the extinction curve, and R(V) has the usual definition of $R(V)\equiv \frac{A_V}{E(B-V)}$. Previous studies of the Galactic extinction law in the NIR indicated that it could be represented as a universal power law with A(λ) ∝ λ^−β, where the exponent β is constant for any line of sight (the usually agreed upon value of 1.6 ≤ −β ≤ 1.8; see Draine 1989; Sneden et al. 1978; Schultz & Wiemer 1975). FM09 showed that β is not only dependent on the line of sight, but that the wavelength range used to probe the extinction actually affects the value of β, and hence the value of R(V) that is determined. See Sect. 3.2 of FM09 for a full discussion. Owing to its success in determining extinction towards a wide variety of regions, its flexibility, and simplicity, we have adopted the form of the extinction curve of FM09 using formula 2. This extinction curve mimics a power law whose exponent increases with wavelength and can vary for different lines of sight, while keeping the number of free parameters to just two–R(V) and α. Here the exponent α is substituted for β simply to highlight that it is no longer constant. By using the FM09 extinction curve, which has been characterised using only broadband photometry, we have assumed that it is appropriate for our wavelength range. This can be confirmed in the coming years with JWST spectroscopy.

4.4 The α-R(V) relation

Photometric extinction studies using high mass stars utilise both an expected colour excess E(B – V) and expected brightness M for each star. Both of these values are used in FM09 to determine a and R(V). For nebular observations, while Case B recombination theory can provide E(B – V) via Eq. (1), there is no expected nebular brightness. This prevents us from directly determining R(V) in Eq. (2). FM09 found that α and R(V) are related, which allows us to determine R(V) indirectly through α. To do this, we took the values and 1σ uncertainties for α and R(V) from FM09, and ran Monte Carlo simulations to fit a line between the data points, and determine the uncertainty of the fitting parameters. The α–R(V) relation is shown in Fig. 4. The uncertainties of R(V) and α are correlated, as indicated by the direction of the error bars in Fig. 4. We took care to only allow the data points to vary along the directions indicated by their error bars. We performed the Monte Carlo simulations multiple times, fitting different types of curves to the data. These included a first-and second-order polynomial, as well as a power law. We found that these choices did not significantly affect the final result. The goodness of fits for each curve were very similar, and the best fitting curve closely resembled a straight line in all three cases. Hence, we chose to continue with the first-order polynomial. The equation of the line of best fit is $$\alpha =-0.287(\pm 0.018)\times R(V)+3.309(\pm 0.07).$$(3)

Using this relation, we replaced α in Eq. (1) with Eq. (3), reducing it to a single parameter equation. Having done this, we were now able to begin fitting different extinction curves to solve for R(V).

Fig. 4 Relationship between R(V) and α determined by Fitzpatrick & Massa (2009). The uncertainties in α and R(V) are correlated, and the direction and strength of the error bars reflect that. The orange shaded band is the 1σ uncertainty of the fit, determined by performing 1000 Monte Carlo iterations.

4.5 Determining R(V) and E(B–V)

Following Eq. (1), in order to calculate E(B – V), the decrement: log(r_obs/r_įnt) is divided by the extinction curve difference: κ(λ_line1) − κ(λ_line2). In order for each line ratio to produce the same value of E(B – V), such that the dispersion is zero, the value of the decrements must be equal to, or be a constant multiple of the extinction curve differences. As such, by varying R(V), the dispersion in E(B – V) changes accordingly. By minimising the dispersion, the best fitting value of R(V) can be found. The process of minimising the dispersion in E(B – V) is analogous to performing the χ² statistic test and minimising χ², where ${\chi }^2\equiv \text{Σ}\frac{{({\text{data}}_{\text{i}}-{\text{model}}_{\text{i}})}^2}{{\text{model}}_{\text{i}}}$.

Using Eq. (2), we produced 100 extinction curves with R(V) going from 2 to 8. To determine the best fitting value of R(V), each decrement was fit with all 100 extinction curves, with each extinction curve being normalised to the level of the first decrement (Br₁₁/Br₆). By normalising, any offsets in the y-axis between the extinction curves are removed, and the difference between the decrements (data) and the extinction curve (model) are set entirely by the value of the exponent α in Eq. (2), which we express in terms of R(V) using Eq. (3). Through this, by varying the value of R(V), we could find the value that produced the minimum χ² statistic.

To determine R(V) for each decrement, we ran a Monte Carlo simulation, allowing the value of the decrements to vary within their uncertainties. In order to fully sample the uncertainties of the decrements, we needed to generate 10 000 realisations per spectrum. This larger number of iterations compared to Sect. 4.2.2 is due to the decrement uncertainties being larger than the line flux uncertainties. With each iteration, all 100 normalised extinction curves were fit to the current realisation of the four decrements, calculating the χ² statistic each time. The minimum χ² and corresponding value of R(V) were saved before moving on to the next realisation. This resulted in 471 × 10000 χ²s. We then calculated the median χ² and R(V) for each of the 471 spectra. For all 471 nebular spectra, a typical value of R(V) = 4.9 ± 1.5 was found.

After careful inspection of the decrements with respect to their best fitting extinction curves, we chose to filter our sample of 471 lines of sight. This was based on a combination of low S/N spectra (emission lines with S/N < 3), resulting in large line ratio uncertainties, and spectra that produced unphysical jumps/spikes in the decrement trend. The reason for these jumps/spikes is still not clear to us, but may originate from unresolved issues with the NIRSpec calibration. We ordered the spectra in terms of increasing χ², and chose to consider the best 200. In general, beyond this point the decrements exhibited more unphysical spikes/jumps, with uncertainties becoming very large, of order ≥50%. An example of “reliable” and “unreliable” decrements are shown in Fig. 5. From the 200 most reliable decrements, we found the typical value of R(V) = 4.8 ± 1.06. The nominal value of R(V) is hardly altered by this filtering procedure, indicating that we have not introduced biases into R(V) by only considering certain spectra. We also explored whether the smallest χ² values corresponded to the brightest nebular spectra (with the highest S/N), as there is likely a correlation between brighter spectra and regions with lower extinction/colour excess. However, we did not find any significant correlation between the S/N of the spectra and their resulting χ². As such, we felt it was justified to exclude the sight lines that increased our uncertainty in R(V).

The χ² statistic was used in this analysis as it is a familiar and readily understood measure of goodness-of-fit. However, the minimisation of dispersion in E(B – V) has a more immediate physical meaning, in that the colour excess along the line of sight should be constant, and hence dispersion should be zero when the appropriate extinction curve is used. We performed the above analysis a second time, minimising the dispersion in E(B – V) rather than χ², and filtered our spectra based on increasing values of dispersion. By minimising the dispersion in E(B – V), a typical value of R(V) = 4.97 ± 1.1 was found, and 198/200 spectra were overlapping between the two approaches, showing that the methods are effectively equivalent.

4.5.1 Determining E(B–V)

To determine the colour excess for each spectrum, we fit our 200 spectra with their best fitting extinction curves. This was done via Monte Carlo as before, to propagate the uncertainties of the decrements and R(V) to E(B – V). This produced four values of E(B – V) per spectrum. The median of these four values was calculated to determine the final colour excess for each spectrum. The median for all spectra was calculated, representing the typical value of E(B – V) for our 200 sight lines. We found a typical value of E(B – V) = 0.64 ± 0.27.

Fig. 5 Decrements are shown as black solid lines, with a reliable fit for R(V) = 5.0 (top) and an unreliable fit for R(V) = 2.8 (bottom). The uncertainties of each decrement are indicated with black error bars. Four normalised extinction curves are shown. R(V) = 2 is represented by the blue, sparsely dotted line. R(V) = 3 is the magenta dot-dashed line. R(V) = 5 is the green dashed line, and R(V) = 8 is the red dotted line. Due to the relatively small differences in slope between extinction curves with R(V) = 2 and R(V) = 8 at these wavelengths, small uncertainties in the line ratios are required in order to reliably determine the best fitting value of R(V).

4.5.2 The total extinction A(V)

Given the definition of $R(V)\equiv \frac{A(V)}{E(B-V)}$, having determined both R(V) and E(B – V), we were able to calculate the total extinction A(V). We have found A(V) = 3.1 ± 1.5 mag. This value represents both the foreground extinction due to the diffuse ISM, as well as the local extinction happening within NGC 3603. In order to separate these two distinct regimes of extinction, we have used the Gaia data release 3 catalogue Brown et al. (2021) of extinction values for 95 sources located between 0.12°–0.5° of NGC 3603, at a distance of 6–8 kpc in order to determine the foreground extinction due to the diffuse ISM. We have obtained this data directly from ESA Sky. The visual extinction in the Gaia system, A₀, has a wavelength of 0.547 µm, just slightly shorter than A(V) at 0.551 µm. Due to this small difference in wavelength we consider A₀ and A(V) directly comparable. We have found a foreground extinction of A₀ = 1.96 mag. No measurement uncertainties for A₀ are provided in the Gaia catalogue. This gives a local extinction value of A(V)_local = 1.1 ± 1.5 mag. This value is consistent within its uncertainties with zero. This suggests that all of the extinction that we have measured towards the region may be occurring in the foreground, and that the cluster itself does not contribute significantly towards the total extinction.

4.6 Stacked nebular spectrum

A high S/N, stacked nebular spectrum has been created¹ by combining and averaging all of the usable nebular spectra together, shown in Fig. B.1. This process made use of all 471 uncontaminated spectra. The spectra were first normalised, setting the continuum level to unity. They were then combined together by averaging each wavelength bin. In cases where the flux of a spectrum in a given wavelength bin deviated by ≥10% of the median value, the wavelength bin for that spectrum was rejected. This removed noisy pixels, cosmic rays, and other spikes present in the spectra. The result was an average nebular spectrum for NGC 3603, with a typical continuum S/N = 70. In this spectrum, the Paschen, Brackett and Pfund series are prominently in emission. The Pfund discontinuity can be seen clearly in this spectrum at λ ≥ 2.3 µm. He I emission lines are also prominent throughout the spectrum, along with some weaker HeII emission. The ions Kr III at λ = 2.1986 µm and Se IV at λ = 2.287 µm, produced by the s-process in Asymptotic Giant Branch stars are strongly detected in the nebular spectra. These ions are commonly seen in planetary nebula, (e.g. Sterling et al. 2017). There are also more than 40 molecular hydrogen H₂ lines present, of which a selection have been annotated in Fig. B.1. These lines are only seen in a subset of the nebular spectra. Typically these spectra correspond to regions of enhanced gas density, such as the south-east pillar shown in Fig. 1, though the lines are weakly detected in more diffuse regions as well.

This spectrum has been presented to showcase the rich NIR emission lines produced by the HII region, and we hope that it will be a useful template for identifying emission lines from other nebulae/HII regions.

5 Discussion

5.1 Discussion of assumptions and uncertainties

In this study we have made several assumptions in order to determine R(V) and E(B – V) using formula 2. Here a discussion of those uncertainties is given.

5.1.1 Should E(B-V) be constant for each line ratio?

The basis of our analysis relies upon the correct extinction curve producing a constant value of E(B – V), regardless of which line ratio is used to determine it. In Bautista et al. (1995; B95), the authors measured the extinction A(V) in Orion, using the nebular hydrogen emission lines from the Balmer, Paschen and Brackett series. They found that the value of A(V) increased systematically with the mean wavelength of the series, with A(V) = 1.03 ± 0.05, A(V) = 1.22 ± 0.10, and A(V) = 1.86 ± 0.18 for Balmer, Paschen and Brackett respectively. The authors suggest that this may be a result of optical effects due to the complex dust cloud geometry of the region. We also speculate that it may be a result of higher order spectral series (eg. Brackett, Pfund, Humphrey) coming preferentially from hotter/denser regions of the nebula, where the extinction properties may change. In (B95), the authors go on to say that the increase in A(V) only becomes significant when the Balmer series is compared to the Brackett series (difference in A(V) > 3σ). For a single value of R(V), the change in A(V) must be accompanied by a proportional change in E(B – V). Based on this, if spectral series separated by large wavelengths are used to probe the extinction properties of a region, then additional care must be taken to derive R(V). We feel it is justified in our case to assume a constant E(B – V) for each line of sight, as we make use of just a single spectral series.

5.1.2 Calzetti equations assumptions

We have used Eq. (2) to calculate E(B – V), which itself relies upon two assumptions. These are: (1) the extinction that we are measuring is a purely foreground, uniform screen of dust, and (2) that our extinction law is appropriate. We know that (1) may not be completely valid in our case. The morphology of the gas and dust in NGC 3603 is visibly clumpy and irregular (though the intra-cluster extinction is consistent with zero from our results, albeit with large uncertainties). For a uniform foreground dust screen, all scattered and absorbed light is removed from the line of sight, with none of the photons being redirected back towards the telescope. The more complex geometry of dust clouds within NGC 3603 region means that scattered light can indeed be redirected back towards the telescope, including of course the emission lines that we use to measure the colour excess. We may therefore be measuring photons that have experienced scattering by dust, with blue photons being preferentially scattered compared to red photons, and obtain an underestimate of the amount of colour excess. The effect of scattering on extinction and colour excess due to different dust cloud geometries, degree of clumpiness, and aperture size has been investigated by Krügel (2009; K09). One of their final conclusions is that scattering has a non-negligible impact on colour excess and extinction when the scattering optical depth τ_sca within the telescope aperture is larger than the relative observational error. The expression for τ_sca is given in K09: ${\tau }_{\text{sca}}=\frac{1}{2}\varphi Dn{\kappa }_{\text{sca}}$ where ϕ is the spatial resolution of the telescope in radians, D is the distance to the region in cm, n is the number density in cm⁻³ and κ_sca is the scattering coefficient. Using this expression we can show that thanks to the exquisite resolution of JWST along with the relative proximity of NGC 3603 (compared to extra-galactic sources), scattering should not significantly impact our results. For JWST at λ = 2 µm, ϕ = 0.064″ = 3.15 × 10⁻⁷ radians. The distance to NGC 3603 is D = 2.16 × 10²² cm, the number density is 10000 cm⁻³ ≤ n ≤ 100 cm⁻³ and the standard Milky Way value of κ_sca ≈ 2.5 × 10⁻²²cm⁻². This gives 8.4 × 10⁻⁵ ≤ τ_sca ≤ 8.4 × 10⁻³. The upper limit is approximately an order of magnitude smaller than our relative observational error of ~6 × 10⁻², where our relative observational error refers to the typical relative uncertainty of each decrement. This indicates that the effect of scattering does not significantly affect our results. We acknowledge that without a complete radiative transfer computation, we cannot precisely quantify the effect of scattered light, but based on the estimate of τ_sca above, the impact should not be significant for JWST observations of NGC 3603.

For assumption (2), rather than adopting an extinction curve from the literature, we have tested a range of extinction curves, and so we are confident that this second assumption is satisfied in our analysis. A final point to note on the determination of E(B – V), is that NGC 3603 has a very low Galactic latitude ~−0.5° (Fukui et al. 2013), lying well within the Galactic plane. As such, it is possible that other Galactic star clusters may be incidentally located behind NGC 3603, which could affect the extinction properties that we determine. However, given that NGC 3603 is the optically brightest HII region in the galaxy, we feel it is a safe assumption that it is the dominant source of line emission along our lines of sight.

5.2 Comparing E(B–V) with the literature

The value of E(B – V) = 0.64 ± 0.27 that we have found is significantly lower than other estimates in the literature. In Mof83, they found E(B – V) = 1 .44 ± 0.09 by spectrally classifying nine sources and using existing photometry to determine the colour excess. In Mel89, they found E(B – V) = 1.44 ± 0.13 within 50″ of the central OB cluster. In that study they determine E(B – V) with photometric data using the pair method, also known as the colour-difference method or MK method. This technique compares the observed colour and the theoretical colour of zero age main sequence stars that have been spectrally classified. They obtained a similar result using their nebular spectral data using the recombination lines, but critically, they assumed a Galactic extinction law with R(V) = 3.2. In Pan00 they found E(B – V) = 1.41 ± 0.13 for sources within 50″ of the central cluster utilising the P-Q method from Ogura & Ishida (1981). They derive a value of R(V) = 4.3 ± 0.2 for NGC 3603. When they included NIR photometry from Eisenhauer et al. (1998) the value of R(V) increased to 4.2 ± 0.7, 4.8 ± 0.8, and 5.2 ± 1.0, for the J, H and K filters respectively. The uncertainty was much higher for these NIR observations, and so Pan00 adopt their own value of R(V) = 4.3 ± 0.2 owing to its lower uncertainty. In Sun04, it was found that E(B – V) = 1.25 using the colour-difference method, and using the relation (4) from Guetter & Vrba (1989), $$R(V)=2.45\times \frac{E(V-I)}{E(B-V)},$$(4)

they found a value of R(V) = 3.55 ± 0.12. They also determined a distance modulus V₀ − M_V = 14.2 ± 0.2 for NGC 3603 based on their estimate of the colour excess. Pang et al. (2016) (Pang16) also found a value of R(V) that is higher than the typical galactic value. Within the filter system for HST, they found R_ƒ555w = 3.75 ± 0.87 and a value of E(ƒ435W − ƒ555W) = 1.33 ± 0.12 mag. They do this using the pair method to determine E(ƒ435W − ƒ555W), and perform a similar step to determine the total extinction A_ƒ555W by comparing the measured magnitude and theoretical magnitude of zero age main sequence stars. This allows them to directly calculate R_ƒ555W from the relation $R_{f555W}\equiv \frac{A_{f555W}}{E(f435W-f555W)}$. We note that they rely upon the distance modulus V₀ − M_V from Sun04 to determine the total extinction, but the distance modulus from Sun04 was determined based on a different colour excess law than that found by (Pang16). See Melena et al. (2008) for a full summary of selective extinction values from literature for NGC 3603.

5.3 Spatial sampling of nebular spectra

A possible explanation behind why our value of E(B – V) is so much lower than what has been found in previous studies may be the spatial distribution of the nebular spectra in our sample. The extinction studies mentioned in Sect. 5.2 make use of the colours of massive stars that are close to the cluster centre (≤50″). Our nebular spectra sample a wider area, with the majority of our spectra coming from regions ≥60″ from the cluster core, and some spectra coming from as far as ~ 100″ from the cluster core. Many of our lines of sight sample denser regions of the nebula, where large clouds of gas and dust can be seen. This is qualitatively different from the inner 50″ of the cluster, where the density of gas and dust is much lower. This is evident in Fig. 1, in the top left of the image, where the nebular emission drops off sharply beyond the south-east pillar. This variation in gas and dust density could lead to different extinction properties across the region. Indeed, we have found that the best fitting value of R(V) and resulting E(B – V) varies across the region, with a 1σ spread of ±1.25 and ±0.6 respectively. As such, we believe it is likely that our lower typical value of E(B – V) is not contradicting previous studies, rather it is reflective of the different regions of the nebula that we are probing.

6 Conclusions

Using the NIR nebular spectra of NGC 3603, we have determined a new extinction law towards the giant HII region, mostly corresponding to regions ≥50″ from the cluster centre. We have done this by exploiting the strong Brackett recombination lines present in our nebular spectra obtained with JWST NIRSpec. By fitting a range of extinction curves to each spectrum, we were able to home in on the most appropriate extinction curve, which has the form of FM09’s NIR extinction law. We have found a typical value of R(V) = 4.8 ± 1.06 for the region as a whole. This corresponds to a colour excess of E(B – V) = 0.64 ± 0.27, which is about half the value published in previous works. Our findings diverge from what has been found previously, but we believe this is a result of the spatial sampling of the region, and highlights how spatially variable the extinction law is across the HII region.

We have also created a stacked nebular spectrum by averaging together all of the usable individual nebular spectra, resulting in a typical continuum S/N = 70. This stacked spectrum highlights the emission lines produced by the nebula of NGC 3603, and can be used as a template for identifying otherwise hard to detect emission lines in fainter nebulae.

Acknowledgements

We thank our referee for their suggestions and comments which have greatly improved both the analysis presented here, as well as the readability and flow of the paper.

Appendix A Impact of contaminant stars

Of the 600 nebular spectra obtained with NIRSpec, 129 were not usable due to contamination by an incident star, leaving us with 471 nebular spectra. As mentioned in section 2, a limitation on the number of sources that could be observed in a single MSA configuration came from the intense crowding of NGC 3603. In order to ensure that only well isolated sources were observed, with clean neighbouring background shutters, a criterion was implemented which specified that a source could be considered for observation so long as any other nearby stars were at least five magnitudes less bright in the V band. This approach was successful in so far that every target star had at least one clean background shutter. However, the presence of contaminant stars, even faint ones, has a significant impact on the extracted nebular background spectrum. Figure A.1 shows an extracted spectrum from a nebular background shutter that contained a contaminant star, along with a spectrum from a clean neighbouring shutter. In figure A.1 the upper blue spectrum is an example of a contaminated nebular spectrum. The lower green spectrum is an example of a clean nebular spectrum, with no contamination. Both of these spectra come from neighbouring micro-shutters, and so they should look quite similar to each other in terms of over all brightness and spectral slope. The most obvious difference between the clean and the contaminated spectra is the shape of the continuum, which for the contaminated spectrum, takes on a quasi-periodic oscillating profile. It is not clear to us why this happens, however it is apparent that this occurs only in cases where there is a contaminant star. This phenomenon was also seen in previous work with simulated NIRSpec data Rogers et al. (2022). With the simulated data, contaminants could be added or removed at will, and the same oscillating pattern disappeared from the nebular spectra once the contaminants were removed. It is also clear in the 2D spectrum of affected nebular spectra, where a narrow, bright spectral trace can be seen overlaid on top of the extended nebular emission.

We have attempted extracting contaminated spectra in different ways to try to minimise or remove the effect of the contaminant. NIPS allows us to extract the spectrum in multiple ways. The traditional approach is to simply sum the flux in each detector column, but it is also possible to take a median value of the flux in each column, and multiply this value by the number of pixels per column (in our case for the 2D rectified spectra, this is seven pixels). The benefit of using the median approach for nebular spectra is that, under the assumption that the nebula is spatially uniform across the micro-shutter, the S/N of the extracted spectrum can be improved. This is because there is some slit loss around the edges of the shutter which can essentially be corrected by taking the median flux value. In the case where the contaminant spectrum is located near the upper or lower edge of the spectral trace, the median pixel value would correspond to the nebular light, excluding the contaminant flux. In the more common case of bright contaminant stars, whose spectral trace can cover ≥ 50% of the seven pixels per column, the stellar flux dominates the spectrum, and taking the median effectively considers only the stellar flux. For this reason we opted to stick with the traditional sum method.

Another approach to remove the contaminant flux is to simply not count the detector rows that contain the stellar flux. Rectifying the spectrum means that the spectral trace runs perfectly parallel with the dispersion direction, making it relatively simple to identify which rows contain the contaminant flux, and skip them. Unfortunately, the resampling of the flux during rectification does not work perfectly, with some pixels receiving more flux than they should, and others receiving too little. The result is that extracting only a portion of the nebular spectrum - the portion that does not contain contaminant flux, results in its own oscillating continuum due to the resampling. This is most easily seen in data like ours, with high S/N and high resolution. Based on this, we opted to simply exclude any nebular spectrum that contained a contaminant. This brought the number of usable nebular spectra from 600 to 471.

Fig. A.1 Two nebular spectra, one that has been extracted from a clean nebular region free of contaminants, and another that is dominated by a contaminant.

Appendix B Stacked nebular spectrum

Fig. B.1 Stacked spectrum, resulting from the averaging and combining of the individual nebular spectra.

References

All Tables

All Figures

Fig. 1 Region of NGC 3603 with the projected micro-shutter’s angle and dimensions. This RGB image is composed of High Acuity Wide-field K-band Imager (HAWK-I) F212n H2 (Red), HST F160W (Green), and HST F168n (Blue). The stellar sources not discussed in this paper lie within the central micro-shutters, while the upper and lower shutters observe the nebula. The zoomed in section containing the source 3796 illustrates the nodding pattern utilised. JWST nods up and down, moving three neighbouring open micro-shutters north-east and south-west, observing 6 nebular regions in total. The position of the OB cluster is indicated by the blue star and arrows.

In the text

	Fig. 2 Example of four nebular spectra exhibiting a raised continuum level around Br10−4 as a result of the NIRSpec calibration. This effect produces unphysical recombination ratios with respect to Case B recombination, and so this line has not been included in the remainder of the analysis.
In the text

	Fig. 3 Five Brackett series line profiles are shown from a representative spectrum, plotted in log-scale. The best fitting Gaussian profile is shown as a red dashed line.
In the text

	Fig. 4 Relationship between R(V) and α determined by Fitzpatrick & Massa (2009). The uncertainties in α and R(V) are correlated, and the direction and strength of the error bars reflect that. The orange shaded band is the 1σ uncertainty of the fit, determined by performing 1000 Monte Carlo iterations.
In the text

Fig. 5 Decrements are shown as black solid lines, with a reliable fit for R(V) = 5.0 (top) and an unreliable fit for R(V) = 2.8 (bottom). The uncertainties of each decrement are indicated with black error bars. Four normalised extinction curves are shown. R(V) = 2 is represented by the blue, sparsely dotted line. R(V) = 3 is the magenta dot-dashed line. R(V) = 5 is the green dashed line, and R(V) = 8 is the red dotted line. Due to the relatively small differences in slope between extinction curves with R(V) = 2 and R(V) = 8 at these wavelengths, small uncertainties in the line ratios are required in order to reliably determine the best fitting value of R(V).

In the text

	Fig. A.1 Two nebular spectra, one that has been extracted from a clean nebular region free of contaminants, and another that is dominated by a contaminant.
In the text

	Fig. B.1 Stacked spectrum, resulting from the averaging and combining of the individual nebular spectra.
In the text