EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 681 (January 2024) A&A, 681 (2024) A20 Full HTML
Open Access
Issue
A&A
Volume 681, January 2024
Article Number A20
Number of page(s) 16
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202348097
Published online 03 January 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

The 30 Doradus region (also known as the Tarantula nebula) in the Large Magellanic Cloud (LMC) is the most luminous and most energetic star-forming region in the Local Group. In its centre lies the OB association NGC 2070, which itself hosts the compact subcluster R 136 at its core, a young (1–2 Myr; Crowther et al. 2016; Bestenlehner et al. 2020) star cluster that contains some of the most massive stars known (Crowther et al. 2010; Brands et al. 2022). 30 Dor is therefore the closest example of a massive extragalactic star burst, powered by a young massive cluster, and it is often considered as a local analogue of the extreme star and star cluster formation conditions in more distant interacting galaxies (e.g. Bastian et al. 2013) or at high redshift (e.g. Bouwens et al. 2021).

Consequently, understanding how star formation proceeds in 30 Dor has been the focus of many observational campaigns across the full wavelength range. X-ray and ultraviolet studies have established the significant contribution of the R 136 cluster to the ionising flux and mechanical feedback in the region (Bestenlehner et al. 2020; Cheng et al. 2021; Crowther et al. 2022), but early ground-based studies already suggested that star formation is ongoing beyond the central cluster (e.g. Walborn 1991; Hyland et al. 1992; Rubio et al. 1992). For example, Walborn & Blades 1987 and Walborn (1991) identified three O-stars that are deeply embedded in nebular knots to the north-east and west of R 136 with spectroscopy, while Rubio et al. (1998) detected several infrared (IR)-bright sources in the nebular filaments surrounding the central cluster with deep ground-based imaging.

Optical imaging with the Hubble Space Telescope (HST) Wide Field Planetary Camera 2 (WFPC2) revealed highly complex nebular structures with compact dust regions and extended pillars (Hunter et al. 1995b; Scowen et al. 1998; Rubio et al. 1998; Walborn et al. 2002). Complementary to ground-based IR studies, Walborn et al. (1999) presented HST Near-Infrared Camera and Multi-Object Spectrometer (NICMOS) imaging of several fields in 30 Doradus, resolving previously identified IR-bright sources into small clusters with multiple components. Brandner et al. (2001) used the same data to identify pre-mainsequence (PMS) stars.

Using HST Wide Field Camera 3 (WFC3) observations of the 30 Dor nebula, De Marchi et al. (2011b, 2017) identified more than 1150 PMS stars in a 40 × 40 pc region. Similarly, an extended distribution of upper main-sequence (UMS) stars was found by the Hubble Tarantula Treasury Program (HTTP; Sabbi et al. 2013, 2016), a multi-band HST imaging survey spanning a much larger region of 30 Dor. Using these data, Ksoll et al. (2018) identified more than 16000 likely PMS stars with a machine-learning approach. These stars are distributed across a 200 × 175 pc region in loose and sometimes filamentary structures.

Using observations in the mid- and far-IR with the Spitzer and Herschel space telescopes, tens of embedded young stellar objects (YSOs) have been identified based on their colours and spectral energy distributions (e.g. Whitney et al. 2008; Gruendl & Chu 2009; Carlson et al. 2012; Walborn et al. 2013). Recently, Nayak et al. 2023 reanalysed data from the Spitzer SAGE (Meixner et al. 2006) and Herschel HERITAGE (Meixner et al. 2013) surveys to identify a large sample of nearly 300 candidate high-mass YSOs in 30 Doradus, some distributed in clumps 45 pc away from R 136. High angular resolution ALMA CO(2–1) observations of the region presented by Wong et al. (2022) also found a complex network of filamentary structures in cold molecular gas, but Nayak et al. (2023) found only ~40% of their YSO candidates to be associated with the CO molecular gas, suggesting that a majority of the cold H2 gas is not traced by CO.

Overall, observational studies of 30 Doradus have revealed it to be a complex region organised in clumps and filaments, in which multiple episodes of star formation have occurred in the past ~30 Myr and are still ongoing in the molecular gas around 30 Dor (Rubio et al. 1992; Walborn et al. 1999; Kalari et al. 2018; Nayak et al. 2023). R 136 was found to host multiple generations of stars as young as 0.5 Myr (Crowther et al. 2016; Bestenlehner et al. 2020), whereas the Hodge 301 cluster, located 45 pc to the north-west, is substantially older, with an age of up to 30 Myr (Lortet & Testor 1991; Grebel & Chu 2000; Cignoni et al. 2016). Using HST, De Marchi et al. (2017) found low-mass PMS stars as old as 50 Myr to be separated spatially from R136, and widely distributed O and B stars found by the VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011) indicate distributed star formation within the past few million years (Schneider et al. 2018). Additionally, the massive YSO candidates identified from Spitzer and Herschel data by Nayak et al. (2023) suggest an increase in the star formation activity in the region in the past few million years.

In this work, we build upon existing studies of the young and old stellar populations in the 30 Doradus region using JWST NIRCam photometry of the region. While Spitzer and Herschel have been instrumental in identifying deeply embedded YSOs that are invisible to the optical cameras of the HST, their limited sensitivity and angular resolution restricted the analysis to short-lived YSOs with masses of ≳ 5 M⊙. On the other hand, HST-based studies in the optical enabled identification of low-mass PMS stars down to ~0.5 M⊙ that have a considerable age spread (De Marchi et al. 2017), but miss the most embedded phases. With its high sensitivity, high angular resolution, and coverage in the infrared, the JWST now allows us study PMS stars of various ages even in highly dust-obscured regions, as shown, for example, with the JWST NIRCam observation of NGC 346 in the Small Magellanic Cloud (SMC; Jones et al. 2023) or NGC 6822 (Lenkić et al. 2023; Nally et al. 2023). In this work, we used publicly available NIRCam data of the central field from Early Release Observations (ERO programme 2729, PI: K. Pontoppidan) as well as two NIRCam fields separated from the centre of the cluster (programme 1226, PI: De Marchi) to distinguish the various stellar populations in 30 Doradus.

The paper is structured as follows. Section 2 describes the data reduction and photometry workflow, and Sect. 3 describes how we corrected for extinction. Our methods of selecting different subpopulations are presented in Sect. 4. Section 5 describes their spatial distributions in the 30 Doradus region. We discuss our findings with respect to the literature in Sect. 6 and conclude in Sect. 7. Throughout this work, we assume a distance to the LMC and 30 Doradus of D = 51.3 kpc, corresponding to a distance module of µ = 18.55 mag (Panagia et al. 1991; Panagia 2005). At this distance, 1′ corresponds to 14.9 pc and 1″to 0.25 pc, respectively.

thumbnail Fig. 1

Greyscale DSS2 red image of 30 Doradus. The red, orange, and blue rectangles show the footprints of the central, northern parallel, and southern parallel fields, respectively. The black circles indicate the location of well-known clusters in 30 Doradus. North is up, and east is to the left. The greyscale image has a size of 17.5′ × 17.5′ (260 × 260 pc). The centre field spans a 7.4′ × 4.4′ region (110 × 66 pc), and the parallel observations consist of two fields each with an area of 2.2′ × 2.2′ (33 × 33 pc), corresponding to the footprint of the NIRCam A and B modules.

2 Data

In this work, we used publicly available JWST NIRCam data from the Early Release Observations programme 02729 as also described in Fahrion & De Marchi (2023), where we used the JWST data in combination with HTTP HST data to derive the extinction law in 30 Dor from 0.3 to 4.7 µm. Here, we complemented the public JWST data of the central region of 30 Doradus (called central field or centre field in the following) with two additional NIRCam pointings that were observed on 8 July 2023 as parallel observations to NIRSpec observations in programme 01226 (PI: De Marchi). We refer to these fields as the northern and southern parallel fields. Table 1 lists basic details of the observations such as the filter combinations and exposure times. Figure 1 shows an overview of the 30 Doradus region with a DSS2 red image and the footprint of the JWST data highlighted. Figure 2 shows the image from the offical release1.

2.1 Data reduction

For this paper, we further refined our photometric catalogues presented in Fahrion & De Marchi (2023). Firstly, we reran the NIRCam data reduction pipeline steps 2 and 3 using JWST_PIPELINE version 1.10.02 with the CDRS context map JWST_1089.PMAP. We aligned the world coordinate systems to the Gaia DR3 catalogues. The main difference to the older data reduction version used in Fahrion & De Marchi (2023) is the treatment of bad-quality pixels that were set to zero in the old version and are now NaN values. This leads to a smoother appearance in the final mosaics, but does not affect the derived photometry.

thumbnail Fig. 2

False-colour images of the centre field from the official image release. Red: F444W, orange: F335M, green: F200W, and blue: F090W. The field shows a region of 7.4′ × 4.4′ (110 × 66 pc). Credit: NASA, ESA, CSA, STScI, Webb ERO Production Team.
Table 1

Overview of NIRCam observations.

2.2 Source detection and photometry

While Fahrion & De Marchi (2023) used only the bright red clump stars, we are now interested in the full stellar populations. For this reason, we repeated the source extraction and photometry. For source detection, we used the combined mosaics with an optimised background subtraction to improve the detection of faint sources in regions with strong nebular emission. As in Fahrion & De Marchi (2023), we used PHOTUTILS MMM-BACKGROUND to estimate and subtract the 2D background. Then, we detected all sources in the background-subtracted image that have two connected pixels 3σ above the root mean squared of the background. For the narrow-band filters such as F187N and F470N, we used a 2σ threshold. Then, this master list from the mosaics was used to obtain forced-aperture photometry on the individual exposures using an aperture radius of 2.5 pixels. The background was determined in an annulus with an inner radius of 4 and an outer radius of 5.5 pixels and was subtracted from the aperture flux. To ensure that the world coordinate system was consistent between mosaics and individual exposures, we used calibrated but not combined TWEAKREGSTEP files. They were corrected for the 1/f striping using the routine of Chris Willot2.

As another advancement, we now also performed PSF photometry using PHOTUTILS with the WebbPSF models version 1.1.0. We used a grid of 25 PSFs per detector. For the aperture photometry, the PSF photometry was applied to the individual exposures using the master catalogue from the respective mosaic, but the centroid positions were fitted freely. We note that this only leads to small adjustments on a sub-pixel level in the centroid positions and does not lead to source confusion, but optimises the fits. For each exposure, the two-dimensional background was derived using a box size of 5 × 5 pixels and a filter size of 3 pixels. This background was then supplied to the PSF-fitting routine. We also experimented with directly creating PSF models from the data. However, due to the strongly varying background, we found that the WEBBPSF models perform better. The corresponding aperture corrections are listed in Appendix A.

For each filter, we first obtained the photometry in each exposure and then combined this into a master table. We calculated the magnitudes as the mean from each individual exposure and the corresponding uncertainties as the standard deviations. Then, the catalogues for each filter were combined into a master catalogue, considering only sources that were detected in at least two JWST filters. The final catalogue of the centre field has 149 589 sources. Out of these, 109 782 (73%) can be matched to the HTTP catalogue in the NIRCam footprint. Of those without a HST match, more than half are redder than F090W – F200W = 1.5 mag.

For the parallel observations, we used the same data reduction and source extraction steps. In the northern parallel fields, 26 610 sources are found in at least two of the four filters, but only 3465 of them have a match to the HTTP because the overlapping footprints are small. In the southern parallel observations, 37 680 sources are detected, 12 636 of which are matched to the HTTP. Although the exposure time is the same for both parallels, the larger number of filters and the choice of wide filters over medium wide filters in the southern parallel observations leads to the larger number of detected sources.

thumbnail Fig. 3

False-colour images of the northern parallel (top) and southern parallel (bottom) fields. The orange arrow indicates the orientation on the sky. The images were constructed using different JWST NIRCam filters in the respective fields. The colours for the northern parallel field are red: F430M, yellow: F405N, blue: F182M, and cyan: F187N. The colours for the southern parallel field are red: F444W, yellow: F405N, blue: F200W, and cyan: F187N. Each large square has a field of view of 2.2′ × 2.2′(33 × 33 pc). The dark space in the middle corresponds to the 44″ gap between the NIRCam A and B modules, while the cross shape arises from the ~5″ gaps between the four individual detectors in the short-wavelength channel. In dithered observations such as the observations of the centre field, these gaps disappear. The location of these fields on the sky is shown in Fig. 1.

2.3 Colour images of the parallel observations

Figure 3 shows false-colour images of the two parallel fields, created by combining all of the available NIRCam filters. Red represents long wavelengths, and blue represents short-wavelength channels.

These images clearly show nebulosity with a high level of detail in both parallel fields, in particular, in the western parts of the observations, which are closer to the centre of 30 Doradus. The eastern parts also show faint nebular features. The north-west corner of the southern parallel field covers part of the SL 639 cluster, for which Britavskiy et al. (2019) reported an age of ∼22 Myr based on VFTS data (see also Evans et al. 2015).

thumbnail Fig. 4

Extinction correction for the centre field. Top left: original F090W – F200W CMD. The cyan line represents a 4 Myr (log(Age) = 6.65) MIST isochrone with metallicity [Fe/H] = −2.5 dex (solid in the upper main sequence). The red parallelogram shows our selection of UMS stars that we used to derive the extinction by placing them back onto the isochrone. Top right: extinction-corrected CMD. Bottom left: normalised histogram of the derived extinction values AF090W. The mean is at 1.31 ± 0.20 mag, and the arrow in the top left panel corresponds to this extinction. Bottom right: map of derived extinction values. The smooth appearance is caused by the smoothing parameter ϵ.

3 Correcting for extinction

We derived the extinction law in 30 Doradus in Fahrion & De Marchi (2023) using the same JWST NIRCam data in the centre field as well as archival HST data from the HTTP. We now used the derived law to correct the colour-magnitude diagrams (CMDs) of the central field for extinction using a similar approach as described in Ksoll et al. (2018) or De Marchi et al. (2011a, 2017). As we are mainly interested in correcting young populations and PMS stars, we based our extinction correction on UMS stars for which the absolute extinction can easily be derived by placing them back onto a theoretical isochrone. We first selected UMS in the F090W – F200W CMD as shown in Fig. 4 and then placed each UMS star back to its theoretical location on an unreddened 4 Myr MIST3 isochrone (Dotter 2016; Choi et al. 2016) with [Fe/H] = −0.35 dex (Choudhury et al. 2021, see also Sect. 4.1). This isochrone already includes also extinction in the Galactic foreground. For this component, we used the extinction law from Cardelli et al. (1989), with RV = 3.1, and AV = 0.22, which is the value that is generally adopted for this line of sight (e.g. Brunet et al. 1975; Isserstedt 1975; Fitzpatrick & Savage 1984). As the extinction law derived in Fahrion & De Marchi (2023) has only considered the slope of the extinction vector rather than the absolute extinction value, we did not attempt to consider the total extinction correction, which would also include foreground extinction (see e.g. Maíz Apellániz et al. 2014). In this way, we calculated the absolute extinction in F090W for all selected UMS stars across the field of view.

We then used the UMS stars to derive the F090W extinction for the full sample of stars by calculating the distance-weighted extinction of the nth nearest UMS stars (see also Ksoll et al. 2018), (1)

with weights wi, (2)

where di is the distance to the nth nearest UMS star, and ϵ is a smoothing parameter in the same units as di. Using these parameters, we find an mean extinction of AF090W = 1.31 ± 0.20 mag, as shown in the bottom left panel of Fig. 4.

After testing, we chose to use the 20 nearest UMS stars with a smoothing of ϵ = 5″ around each star in the catalogue. To test this correction, we also applied it to the UMS stars, but in this case, we first removed the respective UMS star from the sample and then used other nearby UMS stars to derive its mean extinction. We found that these parameters resulted in a narrow UMS sequence that closely matched the theoretical isochrone. The right panel in Fig. 4 shows the extinction-corrected CMD in F090W – F200W. We used the extinction law from Fahrion & De Marchi (2023) to also correct all stars in the other JWST NIRCam bands.

The extinction map in the bottom right panel of Fig. 4 shows the distribution of the derived extinction values, clearly illustrating that the extinction follows the optical appearance of the nebula. Figure 4 shows that the extinction-corrected CMD fits the 4 Myr isochrone much better and shows a more narrow UMS. The still extended red clump shows that the correction is not perfect, especially for old populations. However, this is not expected either because by choosing to use the UMS stars for correction, we aimed to correct the extinction for the young populations. The old populations such as red giant stars are not spatially correlated with the UMS stars and therefore cannot be dereddened using UMS stars as reference. Additionally, we caution that we made the assumption that PMS stars follow the same extinction law as the UMS, even though the intrinsic colours of these populations are different and non-linearities in the extinction can have an effect, especially at large absolute extinctions (e.g. Maíz Apellániz et al. 2014; Maíz Apellániz & Barbá 2018). Following Maíz Apellániz et al. (2020), at absolute extinctions in F090W of ~ 1.3 mag, we might expect to see effects at the level of 0.05 mag, similar to the typical uncertainty of the PMS stars. Consequently, disentangling this effect with our data is challenging, especially given that the typical scatter in extinction values of neighbouring UMS stars are about 0.1 mag.

We did not apply a similar extinction correction in the parallel fields because the filter choices were made to identify sources with emission line excess (see Sect. 4.2). As a consequence, no colour combinations are available in which the upper main sequence can be well separated from the red giant branch. Additionally, only a small subset of stars in the parallels overlap with the HTTP footprint and have available HST colours that could be used to derive the extinction. However, we note that we are still able to identify PMS stars and YSOs based on their excess emissions, as described below.

thumbnail Fig. 5

Extinction-corrected CMDs of the central field with MIST isochrones between 0.5 and 20 Myr. The different panels show different colour combinations. Only stars with photometric uncertainties smaller than 0.5 mag are plotted. The black arrows indicate the direction of extinction.

4 Selecting subpopulations

In the following, we describe our approaches to select subpopulations of different ages using colour–magnitude and colour–colour diagrams. For the central field, we used extinction-corrected data.

4.1 Comparison with isochrones

We compare the extinction-corrected CMDs of the centre field in Fig. 5 with theoretical MIST isochrones with ages of 1, 4, 8, and 20 Myr. The isochrones use [Fe/H] = −0.35 dex (see Choudhury et al. 2021) and a Milky Way foreground extinction of AV = 0.22 mag following the law from Cardelli et al. (1989).

The CMD shows that isochrones with ages even younger than 1 Myr are needed to describe the red and faint populations, but some stars are significantly redder than any of the isochrones. This is particularly clear in the right panel, which shows the F200W – F444W CMD. Here, a population of sources with a significant excess towards longer wavelengths is visible. Additionally, the comparison between these two panels illustrates that two blue JWST filters with sufficient separations in their central wavelengths (e.g. F090W and F200W) are needed to different populations using isochrones. Because filters like this are missing in the parallel observations, comparisons with isochrones cannot be used there to distinguish between populations.

We used the MIST isochrones to make a simple selection based on the expected ages of different populations in 30 Doradus. Firstly, we again selected the UMS stars in the CMD, this time, using a simple rectangle and a magnitude limit of F090W = 18.5 mag. For the faint populations, we only considered stars fainter than F090W = 19.5 mag to avoid confusion with giant stars. In the same way as for the youngest stars, we selected all stars redder than the 0.5 Myr MIST isochrone, as shown in Fig. 6. PMS stars with intermediate ages were selected to lie in the region bounded by the other isochrones, selecting stars with ages between 0.5 and 1, 1 and 2, 2 and 4, 4 and 8, 8 and 20, and stars older than 20 Myr. For the oldest stars, we chose an even stricter brightness cut of F090W = 20 mag. To account for photometric uncertainties, we only made these selection for stars with uncertainties <0.1 mag in both the F090W and F200W filters. We note that this selection reaches down to F090W = 25 mag, mass of 0.5 M⊙ for a 20 Myr isochrone, 0.25 M⊙ for a 4 Myr isochrone, or 0.1 for a 1 Myr isochrone. This shows that low-mass PMS stars are included. At the same time, the brightness cut at 19.5 mag sets a limit to stars less massive than 3 M⊙.

As the focus of this work is to investigate how JWST NIRCam can be used to distinguish different populations, we did not consider a more refined fitting of ages and masses, as is often done with HST data (e.g. De Marchi et al. 2017). However, we note that especially when the extinction is very high, the reddening correction can be a source of uncertainty in the age determination, and in individual cases, it can move stars from one selection to the next.

thumbnail Fig. 6

Selecting populations of different ages in a extinction-corrected CMD. The lines show the same isochrones as in Fig. 5. The coloured dots show the selected stars of different ages. We only selected stars fainter than F090W = 19.5 mag to avoid confusion with older giant stars.

4.2 Populations with excess emission

The large number of filter combinations that also consider multiple narrow-band filters allows us to select and study stellar populations within 30 Doradus that show an excess compared to the main populations. This excess can either be caused by emission from certain emission lines (e.g. the HI Paα, Brα, or H2 lines), or from warm dust that emits at long wavelengths. We selected sources that show excess emission using colour cuts, as shown in Fig. 7. This figure shows CMDs using the colours described below.

We considered excess towards the mid-infrared (called IR excess in the following) using the F200W – F444W colours (or F182M – F430M for the northern parallel field). To select excess due to emission from the Paα line, we used F200W – F187N (or F182M – F187N) colours, and for excess due to emission from molecular hydrogen, called H2 excess, we used F444W – F470N (for the centre field and southern parallels). Additionally, the parallel observations were covered with the Brα narrow-band filter F405N, so that we used F444W – F405N for the southern and F430M – F405N for the northern parallel field to select sources with excess in this emission line. Finally, the southern parallel observations also used the F212N narrow band, which covers the 2.12 µm line of H2. To select excess sources, we used the F200W – F212N colour in this case.

Our colour cuts are listed in Table 2. We required the selected sources to have larger (redder) colours than these thresholds. To account for the photometric uncertainties, we only considered sources with magnitude uncertainties below 0.1 mag. Additionally, we only selected sources that met the colour cuts with at least 5σ significance, where σ is the combined uncertainty of the respective two filters.

The last column of Table 2 reports the resulting number of stars that meet these criteria. They are plotted in different colours in the CMDs of Fig. 7. Especially the centre field, many of the excess sources also include bright stars, whereas the parallel observations farther away from the centre of 30 Doradus tend to show fewer excess sources and fainter sources, illustrating that especially massive star formation is mainly limited to the central region. Nevertheless, there are still sources with IR or emission line excess in the parallel fields, clearly indicating ongoing or recent star formation in the parallel fields. These sources lie as far as ~100 pc from the centre of 30 Doradus.

We note that visual inspection of the sources with an Paα excess higher than two magnitudes revealed that around half of them are extended sources in the F187N filter. These mini-HII regions are frequently not detected in F090W.

A combined comparison of the different populations in the centre field is shown in Fig. 8. The left panel of this figure reveals that the IR-excess sources are redder than the main populations in this CMD, while the H2-excess sources scatter across the full range of colours, but are mainly relatively bright due to the limited depth of the F470N observations. In some red and bright Paα sources, the Paα emission is probably caused by winds around massive stars (possibly Be stars), but there might also be some massive PMS stars in this selection. Stars with faint Paα excess show a bimodal distribution. Some of them are close to the main sequence, while others lie towards redder colours, in the area of PMS stars. This is similar to the results of De Marchi et al. (2017), who used HST Hα narrow-band observations in conjunction with broad-band filters to identify accreting PMS stars. In this way, they identified bona fide PMS stars with locations both on PMS tracks and close to the main sequence. This was interpreted to show that older PMS stars that have almost reached the main sequence can also still show an excess in hydrogen emission, likely caused by ongoing accretion. With NIRSpec spectroscopy of similar targets in NGC 346, this emission was then confirmed. This clearly shows that PMS stars in the low-metallicity environment of the SMC can accrete for several dozen million years (De Marchi et al. 2023).

thumbnail Fig. 7

Selecting sources with excess emission. Each panel shows a different filter combination that was used to select sources (from left to right) with IR, Paα, H2, Brα, and H2(2.12 µm) excess. The selected sources are shown with coloured dots. As a result of the filter choices, additional influence of reddening is very small.
Table 2

Overview of selected populations with excess emission.

thumbnail Fig. 8

F090W – F200W CMD for the centre field. Left: different excess populations are shown by the coloured dots. Right: triangles show sources that show excess in two or more selections.

4.3 Colour-colour selections

While we made the initial selections based on only one colour, the position of excess sources in colour-colour space gives further insights into their nature. Figure 9 shows different colour combinations for the different fields. We note that combining different colours to create these figures significantly reduces the number of sources, especially when one of the long-wavelength channel narrow-band filters is included because these data are not as deep.

The leftmost panels of this figure show the colour-colour combinations of the centre field, revealing a population of stars that have both excess in the IR and Paα emission. These are interesting candidates for accreting still deeply embedded PMS stars. The right panel in Fig. 8 shows that many of them are indeed located where PMS stars are expected, but some brighter stars also have excess in IR and Paα. These could be candidates of more massive PMS stars. We note again the population of stars with very high Paα excess (F200W – F187N > 2), which also have a significant IR excess. As described above, many of them show extended structures in the F187N filters and are discernible as mini-HII regions in the images.

Similarly, there are stars with both an excess in IR and H2 emission in the centre field. Some of them are also Paα-excess sources. De Marchi et al. (2023) showed with NIRSpec spectroscopy that accreting PMS stars might also show H2 lines from their natal cloud. Inspection of the sources in the F090W – F200W CMD (Fig. 8) finds them to be brighter than our isochrone-based PMS stars, likely because of the limited depth of the long wavelength data. These stars are likely more massive YSOs and appear to fall on the 1–2 Myr isochrones shown in Fig. 5 for example. Interestingly, all sources with excess in both Paα and H2 also show IR excess.

Because there are fewer stars and the distance to the centre of 30 Doradus is larger, the colour-colour diagrams for the northern and southern parallel observations are less strongly populated with excess sources. The southern parallel observations include a handful of sources with both Paα and IR excess emission, as well as some with H2 emission, creating an excess in both the F470N and F212N filters. We find that sources with an excess in Paα do not necessarily also show an excess in Brα. This probably is because Brα is intrinsically fainter than Paα and also because the F444W filter contains the line, hindering the characterisation of the continuum.

thumbnail Fig. 9

Colour-colour diagrams used to select sources with excess emission due to warm dust or emission lines. Left panels: selection of sources in the centre field with Paα (cyan), H2 (magenta), and IR (purple) excess. Middle panels: selection of sources in the northern parallel field with Paα (cyan), Brα (peach), and IR (purple) excess. Right panels: selection of sources in the southern parallel field with Paα (cyan), H2 (4.70 µm, magenta), H2 (2.12 µm, pale green), Brα (peach), and IR (purple) excess.

5 Spatial distributions of different populations

In the following, we focus our analysis on the spatial distribution of the different populations across the field of view. We then also compare to the distribution of cold molecular gas.

5.1 Two-dimensional histograms of spatial distributions

Figure 10 shows the spatial distributions of the various selected populations in the centre field as two-dimensional histograms overplotted on the DSS2 red image for visualisation. Each bin in the histograms corresponds to 15″×15″. Figure 11 shows similar histograms for the parallel fields.

The distribution of UMS stars peaks sharply in the centre of R136, the young massive cluster in the centre of the 30 Doradus region, and shows a rather smooth distribution across the field of view. This underlines the findings by the HTTP and the VFTS that young, massive stars are widely distributed (e.g. Sabbi et al. 2013; Schneider et al. 2018), likely indicating a widely distributed star formation within the last million years. The very central bin on top of the centre of R 136 shows a drop in these PMS stars, however, we find that this because the extremely massive stars in the centre of R 136 are saturated in all of the NIRCam filter and are thus lacking from our catalogue.

In contrast, the youngest and reddest selected PMS stars (<0.5 Myr) show an elongated distribution towards the northeast of R 136, reaching deep within the highly extinguished region. Ground- and space-based NIR detections of PMS stars have already established that star formation is ongoing in this region (e.g. Rubio et al. 1998; Walborn et al. 1999). With increasing age of the isochrones used for selection, this elongated structure becomes weaker. Instead, stars with colours consistent with ages between 1 and 4 Myr are again concentrated on R 136, but not as peaked as the UMS stars. Stars older than 4 Myr begin to show an elongated curved structure that is displaced from the centre of R 136 and instead fills the rather dust-free region in the centre of NGC 2070. This structure is also found in the 8–20 Myr age bin, but disappears for the oldest stars (>20 Myr), where the dusty regions instead become visible as regions of lower density. We found this to be caused by a selection effect because this population mainly contains faint stars that have a lower detection rate in these regions due to extinction. An additional magnitude threshold of F090W < 22.5 mag for these stars results in a smooth distribution, indicating that this population belongs to the LMC field population.

The third row of Fig. 10 shows two-dimensional histograms for sources selected based on their excess emission. The IR-excess sources show a very similar elongated distribution as the youngest PMS stars (<0.5 Myr). The Paα-selected stars reach the highest densities in the central region of 30 Doradus and overlap with the regions that also have a high density of young and IR-excess populations. Additionally, Paα-excess sources appear at the edges of the nebular regions at the borders of the central NIRCam field of view. Sources with H2 excess are concentrated towards the centre of 30 Doradus and appear to trace the molecular filaments (see also Fig. 12). The parallel fields also contain these excess sources, even as far as 100 kpc from R 136 (Fig. 11).

In the bottom row of Fig. 10, we show histograms when two or more selections are combined. We find that sources that both have an IR excess and are selected to be redder than the 8 Myr isochrone in the F090W – F200W CMD are clearly elongated towards the north-east, while IR-excess stars with older ages are more uniformly distributed across the field. This highlights that the youngest PMS stars often also have IR excess, indicating that they have very red spectral energy distributions. This is expected for PMS stars that are still deeply embedded in their natal cloud (e.g. Robitaille et al. 2007).

Sources with Paα excess and ages younger than 8 Myr are found to be located on regions with high nebulosity in the central region, while almost all Paα excess stars with older ages are located in the other regions of the cloud. We suggest that this indicates that the radiative feedback from the massive stars in R 136 is efficient in destroying the accretion disks of PMS stars in the central regions, while PMS stars at larger radii are able to accrete longer. This erosion of accretion disks due to photo-evaporation was also suggested by De Marchi et al. (2011a) based on the finding that older accreting PMS stars are displaced from the centre of R 136 based on Hα narrow-band HST observations of a central field in 30 Doradus.

thumbnail Fig. 10

Spatial distributions of different populations of the centre field as two-dimensional histograms with 15″ bin size overlaid on the DSS2 red greyscale image. Top and first row: selections based on the F090W – F200W colour magnitude diagram and comparisons with isochrones as shown by two-dimensional histograms. Third row: histograms show individual stars selected based on their excess. Bottom row: histograms showing stars that fall in two or more selections. A scale bar of 2′ is shown in the first panel.
thumbnail Fig. 11

Spatial distributions of different populations of the northern (top) and southern parallels (bottom) as two-dimensional histograms with 15″ bin size overlaid on the DSS2 red greyscale image (see Fig. 1 for the placement of these fields in the 30 Doradus region). Histograms show individual stars selected based on their excess.
thumbnail Fig. 12

Spatial distribution of different populations in the centre field on top of the F335M mosaic (top left and as background in each panel). The red arrow indicates the location of the central cluster R 136, and the red crosses mark the positions of previously well-studied IR-bright knots (e.g. Walborn et al. 1999). Top right: 12CO moment zero maps from ALMA (data from Wong et al. 2022) as green to yellow contours. The coloured dots show individual stars selected in different populations. The image has a field of view of 7.4′ × 4.4′ (110 × 66 pc).
thumbnail Fig. 13

Fraction of sources in different selections that overlap with the ALMA CO(2–1) map from Wong et al. (2022). The cross shows the fraction of overlap without any intensity constraint, the plusses show the overlap with regions in which the CO intensity exceeds 1.0 K km s−1. The black symbols refer to the full JWST catalogue, and the grey symbols and bands show the expected values for a uniform distribution of stars across the field of view with both JWST and ALMA coverage. These values are higher than for the old stars (>20 Myr) because only a few are detected in the highly extinguished regions.

5.2 Comparison with cold molecular gas

While Fig. 10 compares the spatial distributions of different selected populations against the appearance of 30 Doradus in the optical, we compared the youngest PMS stars and the excess sources in the centre field to the ALMA 12CO(2–1) moment zero map from Wong et al. (2022). As discussed in detail by these authors, the high-resolution CO maps show a high degree of complexity, with fine filaments organised in a larger bow-tie structure that fans outwards from the centre of R 136. By comparing with the NIRCam images, we find that the F335M mosaic matches the CO maps best, likely because this filter traces the emission of polycyclic aromatic hydrocarbon molecules (PAHs) that closely follow the distribution of the molecular gas.

The comparison to the youngest PMS shows that many of them are located in cold molecular gas filaments, especially towards the north. The overlap for the IR-excess sources is slightly less obvious because many sources outside the CO filaments are also found. Additionally, both populations show clumps and substructure, many of which are associated with features in the CO map. In addition, there appears to be a diffusely distributed component to these populations, especially the IR-excess sources, illustrating that recent star formation in 30 Doradus was not only restricted to the densest regions.

Similarly, many stars with Paα excess are found outside and at the edges of the cold molecular regions, but some are also located in the nodes of high (but not highest) flux density. In contrast, some of the H2-excess sources are found in the highest flux density regions, likely indicating embedded regions. Many of the H2-excess sources are also found outside of the cold molecular gas filaments, however.

To determine the overlap between populations and the CO maps, we determined the number of stars that are located in a non-zero entry of the ALMA CO(2–1) map compared to all stars in the footprint of the CO map as seen in Fig. 13. We also show the fraction of stars that overlap with regions in which the CO flux densities ≥1 K km s−1 as a stricter limit.

The highest fractions are reached for stars that are selected to be younger than 0.5 Myr. Here, 65% of these stars are associated with CO, and 53% even fall in the high flux density regions. This underlines the qualitative finding based on their distribution on the CO maps and suggests that these PMS stars are indeed still associated with the cold molecular gas filaments.

The next-highest fraction is found for the stars selected based on the H2 excess as identified by the F470N narrow-band filter. Even though this filter traces emission from excited molecular hydrogen, it appears to still be a viable tracer of the most embedded regions. Similarly, selecting sources with IR excess from broad-band colours results in a high overlap with the CO maps.

For the other selections, we find that the fraction with CO overlap decreases with age to fractions <20% for the oldest stars. While this small fraction shows that the old population is unlikely to be associated with the CO emission, we caution that this small fraction is also driven by the low luminosity of these stars. As the CO emission overlaps with the regions of highest extinctions, we detect fewer of these faint, old stars in these regions, as also shown in Fig. 10. To illustrate this further, we include in Fig. 13 the values for all stars in the JWST catalogue. Of these, 27% fall onto a non-zero entry of the CO map, and 18% lie in a region with ≥1 K km s−1. When we instead assume a completely uniform distribution of stars, 28% would be associated with non-zero CO and 20% with the higher-intensity regions. With a magnitude limit on the old stars (>20 Myr), we can recover these fractions. This indeed suggests a uniform distribution.

Although this simple comparison cannot account for the three-dimensional distribution in 30 Doradus, it already shows that the youngest populations are not all associated with the dense molecular gas as traced by CO. This could either be an observational effect due to the limited depth of the CO observations, or it might indicate a large fraction of CO-dark H2, as also concluded by Nayak et al. (2023). These authors compared the location of massive YSO candidates to the same CO(2–1) map and reported that only 38% were associated with the CO molecular gas.

In this regard, the sources with H2 excess that are not associated with the CO molecular gas are particularly interesting because the excited H2 gas traced by the 4.70 µm line might be directly associated with an accretion disk around the stars and not with the natal molecular cloud, as confirmed with NIRSpec spectroscopy of some PMS stars in NGC 346 (De Marchi et al. 2013). However, spectroscopy of multiple H2 lines is required to determine the excitation mechanism (e.g. Jones et al. 2022).

6 Discussion

We have presented an photometric analysis of JWST NIRCam data in three different fields in the 30 Doradus star-forming region. Using public data of the central field as well as two NIRCam pointings conducted as NIRSpec parallel observations, we derived PSF photometry for more than 200 000 stars detected in at least two NIRCam filters. We used selected different populations by a) comparing with theoretical isochrones and b) applying colour-cuts to select sources with a colour excess due to emission from warm dust or excited emission lines. We then analysed the spatial distributions of different populations.

6.1 Comparing approaches based on optical HST and JWST data

The high sensitivity and spatial resolution of NIRCam allowed us to analyse CMDs with similar techniques and accuracy as is typically used for HST photometry. We found a good match to expectations from theoretical isochrones after dereddening the data of the centre field, and we used these isochrones to select PMS stars in а F090W – F200W CMD with ages between 0.5 and 20 Myr. We note here that we did not fit for the ages, but rather used their location in the reddening-corrected CMD to select stars that would fall in a given age bin. Our selection reaches F090W ~ 25 mag, selecting stars with masses between ~0.1 and 3 M⊙. This mass range is similar to that studied with HST within the HTTP project (Sabbi et al. 2013, 2016; Ksoll et al. 2018). Similar to the results of De Marchi et al. (2011a), Sabbi et al. (2016), and Ksoll et al. (2018), we found that the young PMS stars are centrally concentrated and sometimes organised in loose clumps and filaments, while older stars are displaced from the core of R 136. The distribution of UMS stars is strongly peaked at the location of R 136 but also shows a widely distributed component. This agrees with the distribution of UMS stars in the HTTP (Sabbi et al. 2013) and the finding of widely distributed О and В stars from the VFTS (Schneider et al. 2018).

As an alternative to using broad-band colours alone, De Marchi et al. (2010) established that low-mass PMS stars can be identified with high confidence by using optical HST colours in combination with Hα narrow-band observations to search for Hα excess emission. Created by the gravitational energy released from infalling matter in an accretion disk around PMS stars, this feature can be clearly associated with low-mass star formation and has now long been used to identify PMS stars in various star-forming regions in the Local Group (e.g. De Marchi et al. 2010, 2013; Tsilia et al. 2023; Vlasblom & De Marchi 2023). Using HST WFC3 data, De Marchi et al. (2011b) and De Marchi et al. (2017) used this approach to study the PMS stars with Hα excess in 30 Doradus, finding them to be centrally concentrated. With the JWST, we could now perform a similar approach by using the excess in Paα as traced by the F187N filter. As also found by De Marchi et al. (2017), the Paα sources appear to have a broad range in ages. Some of them are clearly located in the PMS region of the CMD, while others are close to the main sequence. This agrees with the interpretation that PMS stars older than 20 Myr can still be accretion material. Based on NIRSpec spectroscopy, emission lines including Paα were indeed recently identified in PMS stars in NGC 346 in the SMC (De Marchi et al. 2023), and a similar study is planned for a sample of about 100 PMS stars in 30 Dor.

Using filters also at longer wavelengths, we were further able to select stars from the same data with an excess in infrared colours, thereby tracing stars with excess emission due to warm dust. This allowed us to study these low-mass PMS stars even in the highly dust obscured regions of 30 Doradus. This is better than what was previously possible in the optical with HST. Comparing the distributions of CMD-selected PMS stars and IR-excess sources, we found the latter to be distributed farther out in the optically obscured regions of the nebula. Using NIRCam observations of the SMC star-forming region NGC 346, Jones et al. (2023) also found that the JWST-selected PMS stars closely follow the dust filaments. Compared to this region, the PMS stars in 30 Dor appear to be both more widely distributed and numerous. This might be expected because 30 Doradus is a more massive region and hosts the most energetic starburst event in the Local Group.

thumbnail Fig. 14

Zoom-ins around the three IR-bright knots 1, 2, and 3. Top row: colour images of 5″ × 5″regions. Red: F444W, green: F200W, blue: F090W. Bottom row: greyscale F200W images with selected sources as coloured dots. The larger purple dots refer to IR-excess sources, and the smaller maroon, red, orange, and yellow dots are stars selected with isochrones with ages of 0.5, 1, 2, and 4 Myr, respectively.

6.2 Comparison to NIR studies

Over the past decades, 30 Doradus has been studied extensively in the near-infrared regime using both ground-based telescopes and HST NICMOS. Numerous IR-bright sources have been detected, in particular towards the north-east and in the western filament next to R 136 (e.g. Hyland et al. 1992; Rubio et al. 1992, 1998; Walborn et al. 1999; Brandner et al. 2001). In these regions, we now also identify a large number of both IR excess and young PMS stars. A particular focus was placed in the previous studies on some of the brightest sources: three early O-type stars embedded in nebular knots, called knots 1, 2, and 3 (Walborn & Blades 1987; Walborn 1991). HST NICMOS observations of these regions revealed them to consist of multiple stellar components embedded in complex nebular and dusty structures, and they are surrounded by fainter stars (Walborn et al. 1999; Brandner et al. 2001). Follow-up observations with HST WFPC2 further showed the complex nebular structure and embedded populations around these knots (Walborn et al. 2002).

Figure 14 shows the JWST view of these three regions, presented in small cut-outs around each knot (see Fig. 12 for their position in the field of view). The panels in the bottom row show the position of sources selected with isochrones (up to 4 Myr) and with IR excess. While the originally studied O-stars are saturated in at least one of the JWST filters, many of the fainter bluer sources seen close to them are young PMS stars. The sources that appear orange or red in the RGB images are predominantly IR-excess sources in our selection. Additionally, the long-wavelength filter F444W shown in red in these images reveals several sources that are not visible at shorter wavelengths.

6.3 Comparison to Spitzer and Herschel sources

Previous studies in the infrared with Spitzer or Herschel were limited to bright, massive YSOs, but used similar colour-based approaches as in this work, or a spectral energy distribution fitting to select PMS stars in the embedded stages of star formation (e.g. Whitney et al. 2008; Gruendl & Chu 2009; Carlson et al. 2012). To understand how these studies compare to our results, we matched our photometric catalogues to the 299 YSO candidates from Nayak et al. (2023). This matching is not straightforward because our catalogues contain many sources within the photometric uncertainties of these candidates. Therefore, we first selected all stars within 1″ from the YSO candidates and then performed a match when a source with an F335M magnitude within 1 mag of the IRAC 3.6 µm or with an F444W magnitude within 1 mag of the IRAC 4.5 µm magnitude was reported in Nayak et al. (2023). This resulted in 72 matches out of the 488 sources listed in Nayak et al. (2023), 30 of which were classified as YSOs.

We note here that not all sources of Nayak et al. (2023) overlap the NIRCam footprint, but the main reason for the low number of matches probably is that most of the YSO candidates are so bright that they are saturated in the NIRCam imaging. The matched sources exclusively have 2MASS Ks magnitudes >15, while the brighter sources are not recovered in the NIRCam images. Consequently, in our selections, even the matched YSOs are mostly missed due to their brightness, showing that an analysis with JWST like this is highly complementary to previous studies of higher mass YSOs. Additionally, in some cases, the YSO candidates could not be matched in magnitude even when measurements were available. The reason probably is that in some cases, the JWST imaging reveals that the YSO candidates consist of more than one source (see also e.g. Jones et al. 2022).

6.4 30 Doradus in the context of star cluster formation

As some of the oldest stellar systems in the Universe, ancient globular clusters (GCs) are often used as tracers of galaxy evolution. These tracers allow us to constrain the assembly of galaxies (e.g. Brodie & Strader 2006; Beasley 2020). However, because they are so old, it is unclear so far how they formed. In general, two pathways are discussed, which either propose that GCs formed under conditions unique to the high-redshift Universe, for example in dark matter minihaloes, or that they are simply the product of violent star formation in high-density conditions that are rare in the local Universe (e.g. Forbes et al. 2018). In the latter scenario, young massive clusters (YMCs) have been proposed as local analogues of GCs at high redshifts and seem to have similar properties as GCs in terms of masses and densities (e.g. Portegies Zwart et al. 2010; Adamo & Bastian 2018; Adamo et al. 2020b). However, their formation is rare in the Local Group, and they are found in high numbers in interacting galaxies (e.g. Whitmore et al. 2010; Adamo et al. 2020a).

In this context, NGC 2070 in 30 Doradus is a unique region in which to study the formation of massive star clusters in great detail. While there are other YMCs in the Milky Way or the Small Magellanic Cloud (e.g. Westerlund 1 & 2), NGC 2070 with its subcluster R 136 is by far the most massive YMC. Mass estimates of the full cluster range from ~7 ×104 M⊙ to 5 × 105 M⊙, while the R 136 subcluster has an estimated mass between 2.2 × 104 and 1 × 105 M⊙ (Hunter et al. 1995a; Andersen et al. 2009; Bosch et al. 2009; Cignoni et al. 2015).

R 136 is known to have an age of 0.5–2 Myr based on its CMD, and the older more diffuse population in NGC 2070 has an age of ~3–7 Myr (Selman et al. 1999; Sabbi et al. 2012; Cignoni et al. 2015; Bestenlehner et al. 2020; Brands et al. 2022). Figure 15 shows that where we selected stars within a 10″ distance (2.5 pc, the typical effective radius of a GC) from the centre of R 136, we can reproduce this finding with the JWST photometry. Domínguez et al. (2023) were able to explain the presence of the younger R 136 subcluster within NGC 2070 with N-body simulations, finding that the original feedback from the stars that formed in the first generation was likely not sufficient to fully expel the gas, allowing a second generation to form that was more centrally concentrated. They interpreted this within the frequently discussed multiple population problem of GCs (e.g. Bastian & Lardo 2018).

The elongated structure of very young (<0.5 Myr), IR-bright stars towards the north-east of the cluster might present additional evidence for a still ongoing hierarchical formation of this massive star cluster. As hydrodynamical simulations suggest, small subclusters form first in filaments during the collapse of giant molecular clouds and then grow through hierarchical mergers while still accreting gas to form a more massive star cluster on a timescale of a few million years (e.g. Howard et al. 2018; Lahén et al. 2020; Guszejnov et al. 2022). In this context, the observed CO filaments might correspond to inflowing material that supplies cold molecular gas and the PMS stars embedded within to the central region of the star cluster to sustain star formation. At some point, feedback from massive stars is thought to heat and dispel the molecular gas, slowing down star formation (e.g. Krumholz et al. 2019; Krause et al. 2020). Given the radiation pressure from the young massive stars in R 136, the filamentary CO structure could also be interpreted as stemming from outflows (Wong et al. 2022). Moreover, our finding that the older Paα sources are mainly found at larger distances from R 136 might indicate that the feedback from the massive stars in R 136 already impedes the accretion of PMS stars in the centre. When the first supernovae explode in NGC 2070, the expelled gas and star formation will likely be halted, leaving a mostly gasfree cluster with stars that formed within a very short time span of a few million years.

thumbnail Fig. 15

Colour-magnitude diagram of R 136 (orange dots) on top of the full region (black). The orange dots refer to all stars within 10″ from the centre of R 136. The coloured lines show theoretical MIST isochrones.

7 Conclusions

We presented a JWST NIRCam analysis of the young and old populations in the 30 Doradus star-forming region. We focused on how JWST filters can be used to select different populations, and we analysed their spatial distributions in the region. Our main results are summarised below:

  • We presented PSF photometry for more than 200 000 sources with at least two NIRCam magnitudes, about 150000 in the centre field and the remaining sources in two parallel fields. Including a short-wavelength filter bluer than 2 µm is crucial for separating different populations in CMDs using classical approaches based on isochrones. We found that the same methods that are long established with HST photometry can be applied. The PSF photometry closely matches the expectations from theoretical isochrones after a reddening correction is applied;

  • We used different approaches to separate various populations, such as a comparison with isochrones of different ages and colour-cuts to separate populations that show excess emission due to either warm dust as traced by a long-wavelength broad-band filter or due to line emission from excited neutral and molecular hydrogen as traced by a corresponding narrow-band filter;

  • Using isochrones on top of the F090W – F200W CMD, we selected low-mass sources whose colours are consistent with ages between 0.5 and 20 Myr down to masses ~0.1 M⊙. We found that the youngest sources (<0.5 Myr) show an elongated structure towards the north-east of the region that extends deep into the nebular regions. Stars with ages between 1 and 4 Myr are centrally concentrated on R 136. In contrast, PMS stars with intermediate ages between 4 Myr and 20 Myr were found to be spatially separated from the very centre of R 136, in an elongated structure, and stars older than 20 Myr appear to be uniformly distributed across the 30 Doradus region, indicating that they belong to the LMC field population;

  • Selecting stars that show an excess towards long wavelengths using the F200W – F444W colour, we found a very similar distribution as for the isochrone-selected <0.5 Myr PMS stars in a more elongated structure. This suggests that there is a population of still embedded very young PMS stars to the north-east of the central cluster that possibly still feeds its assembly through hierarchical mergers of small subclusters;

  • Using the narrow-band filter F187N, we selected sources with an excess in Paα as a tracer of ongoing accretion onto the PMS star. These sources are concentrated towards the centre of the cluster, but are also found further outwards in the optically visible nebular features of 30 Doradus. Placing the Paα sources in a colour-magnitude diagram showed that they are not only in the PMS locations, but also closer to the main sequence, revealing that the accretion can persist many tens of million years. The F470N narrow-band filter as an indicator for excess due to H2 emission predominantly selects more massive stars because the depth is limited. These sources appear to be distributed along the dust regions, but are also found to be distributed more widely;

  • While most of the selected PMS and YSOs are found in the central field, we also find sources with excess due to warm dust and emission lines in the parallel fields. This suggests that star formation is ongoing or very recent as far as 100 pc from the centre. The nebulosity seen in the colour images supports this;

  • By comparing with available CO(2–1) maps from ALMA (Wong et al. 2022), we found that most of the youngest PMS stars (<0.5 Myr) are associated with the CO molecular gas, whereas older or IR-excess stars do not always overlap with the CO maps. About half of the H2-excess sources do not appear to be associated with CO emission. We interpreted this as further evidence for a large amount of CO-dark H2 gas.

JWST NIRCam, with its many filter choices and enhanced sensibility, has allowed us to combine methods that were previously used separately to perform a comprehensive study of the star formation in 30 Doradus.

Acknowledgements

We thank an anonymous referee for comments and suggestions that have helped to polish this paper. K.F. acknowledges support through the ESA research fellowship programme. Based on observations with the NASA/ESA Hubble Space Telescope and the NASA/ESA/CSA James Webb Space Telescope, which are operated by AURA, Inc., under NASA contracts NAS5-26555 and NAS 5-03127. This work made use of Astropy (http://www.astropy.org/): a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

Appendix A Comparing aperture and PSF photometry

We list in Table A.1 the derived aperture corrections for our default choice of an 2.5 aperture with a background annulus between 4 and 5.5 pixels, as well as other choices. These corrections were determined by calculating the offset between aperture photometry and PSF photometry for all unsaturated stars brighter than 20 mag in the respective filter.

Table A.1

Aperture corrections in magnitudes for different choices of aperture radius.

References

  1. Adamo, A., & Bastian, N. 2018, in Astrophys. Space Sci. Lib., 424, The Birth of Star Clusters, ed. S. Stahler, 91 [NASA ADS] [CrossRef] [Google Scholar]
  2. Adamo, A., Hollyhead, K., Messa, M., et al. 2020a, MNRAS, 499, 3267 [Google Scholar]
  3. Adamo, A., Zeidler, P., Kruijssen, J. M. D., et al. 2020b, Space Sci. Rev., 216, 69 [Google Scholar]
  4. Andersen, M., Zinnecker, H., Moneti, A., et al. 2009, ApJ, 707, 1347 [Google Scholar]
  5. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A & A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  7. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bastian, N., & Lardo, C. 2018, ARA & A, 56, 83 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bastian, N., Schweizer, F., Goudfrooij, P., Larsen, S. S., & Kissler-Patig, M. 2013, MNRAS, 431, 1252 [NASA ADS] [CrossRef] [Google Scholar]
  10. Beasley, M. A. 2020, in Reviews in Frontiers of Modern Astrophysics; From Space Debris to Cosmology, 245 [CrossRef] [Google Scholar]
  11. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  12. Bosch, G., Terlevich, E., & Terlevich, R. 2009, AJ, 137, 3437 [Google Scholar]
  13. Bouwens, R. J., Illingworth, G. D., van Dokkum, P. G., et al. 2021, AJ, 162, 255 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brandner, W., Grebel, E. K., Barbá, R. H., Walborn, N. R., & Moneti, A. 2001, AJ, 122, 858 [NASA ADS] [CrossRef] [Google Scholar]
  15. Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A & A, 663, A36 [CrossRef] [EDP Sciences] [Google Scholar]
  16. Britavskiy, N., Lennon, D. J., Patrick, L. R., et al. 2019, A & A, 624, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Brodie, J. P., & Strader, J. 2006, ARA & A, 44, 193 [NASA ADS] [CrossRef] [Google Scholar]
  18. Brunet, J. P., Imbert, M., Martin, N., et al. 1975, A & AS, 21, 109 [NASA ADS] [Google Scholar]
  19. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [Google Scholar]
  20. Carlson, L. R., Sewiło, M., Meixner, M., Romita, K. A., & Lawton, B. 2012, A & A, 542, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Cheng, Y., Wang, Q. D., & Lim, S. 2021, MNRAS, 504, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  22. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  23. Choudhury, S., de Grijs, R., Bekki, K., et al. 2021, MNRAS, 507, 4752 [CrossRef] [Google Scholar]
  24. Cignoni, M., Sabbi, E., van der Marel, R. P., et al. 2015, ApJ, 811, 76 [Google Scholar]
  25. Cignoni, M., Sabbi, E., van der Marel, R. P., et al. 2016, ApJ, 833, 154 [NASA ADS] [CrossRef] [Google Scholar]
  26. Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [Google Scholar]
  27. Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624 [Google Scholar]
  28. De Marchi, G., Panagia, N., & Romaniello, M. 2010, ApJ, 715, 1 [NASA ADS] [CrossRef] [Google Scholar]
  29. De Marchi, G., Panagia, N., & Sabbi, E. 2011a, ApJ, 740, 10 [NASA ADS] [CrossRef] [Google Scholar]
  30. De Marchi, G., Paresce, F., Panagia, N., et al. 2011b, ApJ, 739, 27 [Google Scholar]
  31. De Marchi, G., Beccari, G., & Panagia, N. 2013, ApJ, 775, 68 [NASA ADS] [CrossRef] [Google Scholar]
  32. De Marchi, G., Panagia, N., & Beccari, G. 2017, ApJ, 846, 110 [NASA ADS] [CrossRef] [Google Scholar]
  33. Crowther, P. A., Broos, P. S., Townsley, L. K., et al. 2022, MNRAS, 515, 4130 [NASA ADS] [CrossRef] [Google Scholar]
  34. De Marchi, G., Giardino, G., Biazzo, K., et al. 2023, Nat. Astron., submitted [Google Scholar]
  35. Domínguez, R., Pellegrini, E. W., Klessen, R. S., & Rahner, D. 2023, MNRAS, 520, 5600 [CrossRef] [Google Scholar]
  36. Dotter, A. 2016, ApJS, 222, 8 [Google Scholar]
  37. Evans, C. J., Taylor, W. D., Hénault-Brunet, V., et al. 2011, A & A, 530, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Evans, C. J., Kennedy, M. B., Dufton, P. L., et al. 2015, A & A, 574, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Fahrion, K., & De Marchi, G. 2023, A & A, 671, A14 [Google Scholar]
  40. Fitzpatrick, E. L., & Savage, B. D. 1984, ApJ, 279, 578 [NASA ADS] [CrossRef] [Google Scholar]
  41. Forbes, D. A., Bastian, N., Gieles, M., et al. 2018, Proc. Roy. Soc. Lond. Ser. A, 474, 20170616 [NASA ADS] [Google Scholar]
  42. Grebel, E. K., & Chu, Y.-H. 2000, AJ, 119, 787 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gruendl, R. A., & Chu, Y.-H. 2009, ApJS, 184, 172 [Google Scholar]
  44. Guszejnov, D., Markey, C., Offner, S. S. R., et al. 2022, MNRAS, 515, 167 [NASA ADS] [CrossRef] [Google Scholar]
  45. Howard, C. S., Pudritz, R. E., & Harris, W. E. 2018, Nat. Astron., 2, 725 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hunter, D. A., Shaya, E. J., Holtzman, J. A., et al. 1995a, ApJ, 448, 179 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hunter, D. A., Shaya, E. J., Scowen, P., et al. 1995b, ApJ, 444, 758 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hyland, A. R., Straw, S., Jones, T. J., & Gatley, I. 1992, MNRAS, 257, 391 [NASA ADS] [Google Scholar]
  49. Isserstedt, J. 1975, A & A, 41, 175 [NASA ADS] [Google Scholar]
  50. Jones, O. C., Reiter, M., Sanchez-Janssen, R., et al. 2022, MNRAS, 517, 1518 [NASA ADS] [CrossRef] [Google Scholar]
  51. Jones, O. C., Nally, C., Habel, N., et al. 2023, Nat. Astron., 7, 694 [NASA ADS] [CrossRef] [Google Scholar]
  52. Kalari, V. M., Rubio, M., Elmegreen, B. G., et al. 2018, ApJ, 852, 71 [NASA ADS] [CrossRef] [Google Scholar]
  53. Krause, M. G. H., Offner, S. S. R., Charbonnel, C., et al. 2020, Space Sci. Rev., 216, 64 [CrossRef] [Google Scholar]
  54. Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA & A, 57, 227 [NASA ADS] [CrossRef] [Google Scholar]
  55. Ksoll, V. F., Gouliermis, D. A., Klessen, R. S., et al. 2018, MNRAS, 479, 2389 [NASA ADS] [Google Scholar]
  56. Lahén, N., Naab, T., Johansson, P. H., et al. 2020, ApJ, 891, 2 [CrossRef] [Google Scholar]
  57. Lenkić, L., Nally, C., Jones, O. C., et al. 2023, ApJ, submitted [arXiv:2307.15704] [Google Scholar]
  58. Lortet, M. C., & Testor, G. 1991, A & AS, 89, 185 [NASA ADS] [Google Scholar]
  59. Maíz Apellániz, J., & Barbá, R. H. 2018, A & A, 613, A9 [CrossRef] [EDP Sciences] [Google Scholar]
  60. Maíz Apellániz, J., Evans, C. J., Barbá, R. H., et al. 2014, A & A, 564, A63 [CrossRef] [EDP Sciences] [Google Scholar]
  61. Maíz Apellániz, J., Pantaleoni González, M., Barbá, R. H., García-Lario, P., & Nogueras-Lara, F. 2020, MNRAS, 496, 4951 [Google Scholar]
  62. Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 2268 [NASA ADS] [CrossRef] [Google Scholar]
  63. Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 62 [NASA ADS] [CrossRef] [Google Scholar]
  64. Nally, C., Jones, O. C., Lenkić, L., et al. 2023, arXiv e-prints [arXiv:2309.13521] [Google Scholar]
  65. Nayak, O., Green, A., Hirschauer, A. S., et al. 2023, ApJ, 944, 26 [NASA ADS] [CrossRef] [Google Scholar]
  66. Panagia, N. 2005, in IAU Colloq. 192: Cosmic Explosions, On the 10th Anniversary of SN1993J, eds. J.-M. Marcaide, & K. W. Weiler, 99, 585 [NASA ADS] [CrossRef] [Google Scholar]
  67. Panagia, N., Gilmozzi, R., Macchetto, F., Adorf, H. M., & Kirshner, R. P. 1991, ApJ, 380, L23 [NASA ADS] [CrossRef] [Google Scholar]
  68. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA & A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  69. Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328 [Google Scholar]
  70. Rubio, M., Roth, M., & Garcia, J. 1992, A & A, 261, L29 [NASA ADS] [Google Scholar]
  71. Rubio, M., Barbá, R. H., Walborn, N. R., et al. 1998, AJ, 116, 1708 [NASA ADS] [CrossRef] [Google Scholar]
  72. Sabbi, E., Lennon, D. J., Gieles, M., et al. 2012, ApJ, 754, L37 [Google Scholar]
  73. Sabbi, E., Anderson, J., Lennon, D. J., et al. 2013, AJ, 146, 53 [NASA ADS] [CrossRef] [Google Scholar]
  74. Sabbi, E., Lennon, D. J., Anderson, J., et al. 2016, ApJS, 222, 11 [NASA ADS] [CrossRef] [Google Scholar]
  75. Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018, A & A, 618, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Scowen, P. A., Hester, J. J., Sankrit, R., et al. 1998, AJ, 116, 163 [NASA ADS] [CrossRef] [Google Scholar]
  77. Selman, F., Melnick, J., Bosch, G., & Terlevich, R. 1999, A & A, 347, 532 [NASA ADS] [Google Scholar]
  78. Tsilia, S., De Marchi, G., & Panagia, N. 2023, A & A, 675, A203 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Vlasblom, M., & De Marchi, G. 2023, A & A, 675, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Walborn, N. R. 1991, in The Magellanic Clouds, 148, eds. R. Haynes, & D. Milne, 145 [CrossRef] [Google Scholar]
  81. Walborn, N. R., & Blades, J. C. 1987, ApJ, 323, L65 [NASA ADS] [CrossRef] [Google Scholar]
  82. Walborn, N. R., Barbá, R. H., Brandner, W., et al. 1999, AJ, 117, 225 [Google Scholar]
  83. Walborn, N. R., Maíz-Apellániz, J., & Barbá, R. H. 2002, AJ, 124, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  84. Walborn, N. R., Barbá, R. H., & Sewiło, M. M. 2013, AJ, 145, 98 [Google Scholar]
  85. Whitmore, B. C., Chandar, R., Schweizer, F., et al. 2010, AJ, 140, 75 [NASA ADS] [CrossRef] [Google Scholar]
  86. Whitney, B. A., Sewilo, M., Indebetouw, R., et al. 2008, AJ, 136, 18 [NASA ADS] [CrossRef] [Google Scholar]
  87. Wong, T., Oudshoorn, L., Sofovich, E., et al. 2022, ApJ, 932, 47 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Overview of NIRCam observations.

In the text
Table 2

Overview of selected populations with excess emission.

In the text
Table A.1

Aperture corrections in magnitudes for different choices of aperture radius.

In the text

All Figures

thumbnail Fig. 1

Greyscale DSS2 red image of 30 Doradus. The red, orange, and blue rectangles show the footprints of the central, northern parallel, and southern parallel fields, respectively. The black circles indicate the location of well-known clusters in 30 Doradus. North is up, and east is to the left. The greyscale image has a size of 17.5′ × 17.5′ (260 × 260 pc). The centre field spans a 7.4′ × 4.4′ region (110 × 66 pc), and the parallel observations consist of two fields each with an area of 2.2′ × 2.2′ (33 × 33 pc), corresponding to the footprint of the NIRCam A and B modules.
In the text
thumbnail Fig. 2

False-colour images of the centre field from the official image release. Red: F444W, orange: F335M, green: F200W, and blue: F090W. The field shows a region of 7.4′ × 4.4′ (110 × 66 pc). Credit: NASA, ESA, CSA, STScI, Webb ERO Production Team.
In the text
thumbnail Fig. 3

False-colour images of the northern parallel (top) and southern parallel (bottom) fields. The orange arrow indicates the orientation on the sky. The images were constructed using different JWST NIRCam filters in the respective fields. The colours for the northern parallel field are red: F430M, yellow: F405N, blue: F182M, and cyan: F187N. The colours for the southern parallel field are red: F444W, yellow: F405N, blue: F200W, and cyan: F187N. Each large square has a field of view of 2.2′ × 2.2′(33 × 33 pc). The dark space in the middle corresponds to the 44″ gap between the NIRCam A and B modules, while the cross shape arises from the ~5″ gaps between the four individual detectors in the short-wavelength channel. In dithered observations such as the observations of the centre field, these gaps disappear. The location of these fields on the sky is shown in Fig. 1.
In the text
thumbnail Fig. 4

Extinction correction for the centre field. Top left: original F090W – F200W CMD. The cyan line represents a 4 Myr (log(Age) = 6.65) MIST isochrone with metallicity [Fe/H] = −2.5 dex (solid in the upper main sequence). The red parallelogram shows our selection of UMS stars that we used to derive the extinction by placing them back onto the isochrone. Top right: extinction-corrected CMD. Bottom left: normalised histogram of the derived extinction values AF090W. The mean is at 1.31 ± 0.20 mag, and the arrow in the top left panel corresponds to this extinction. Bottom right: map of derived extinction values. The smooth appearance is caused by the smoothing parameter ϵ.
In the text
thumbnail Fig. 5

Extinction-corrected CMDs of the central field with MIST isochrones between 0.5 and 20 Myr. The different panels show different colour combinations. Only stars with photometric uncertainties smaller than 0.5 mag are plotted. The black arrows indicate the direction of extinction.
In the text
thumbnail Fig. 6

Selecting populations of different ages in a extinction-corrected CMD. The lines show the same isochrones as in Fig. 5. The coloured dots show the selected stars of different ages. We only selected stars fainter than F090W = 19.5 mag to avoid confusion with older giant stars.
In the text
thumbnail Fig. 7

Selecting sources with excess emission. Each panel shows a different filter combination that was used to select sources (from left to right) with IR, Paα, H2, Brα, and H2(2.12 µm) excess. The selected sources are shown with coloured dots. As a result of the filter choices, additional influence of reddening is very small.
In the text
thumbnail Fig. 8

F090W – F200W CMD for the centre field. Left: different excess populations are shown by the coloured dots. Right: triangles show sources that show excess in two or more selections.
In the text
thumbnail Fig. 9

Colour-colour diagrams used to select sources with excess emission due to warm dust or emission lines. Left panels: selection of sources in the centre field with Paα (cyan), H2 (magenta), and IR (purple) excess. Middle panels: selection of sources in the northern parallel field with Paα (cyan), Brα (peach), and IR (purple) excess. Right panels: selection of sources in the southern parallel field with Paα (cyan), H2 (4.70 µm, magenta), H2 (2.12 µm, pale green), Brα (peach), and IR (purple) excess.
In the text
thumbnail Fig. 10

Spatial distributions of different populations of the centre field as two-dimensional histograms with 15″ bin size overlaid on the DSS2 red greyscale image. Top and first row: selections based on the F090W – F200W colour magnitude diagram and comparisons with isochrones as shown by two-dimensional histograms. Third row: histograms show individual stars selected based on their excess. Bottom row: histograms showing stars that fall in two or more selections. A scale bar of 2′ is shown in the first panel.
In the text
thumbnail Fig. 11

Spatial distributions of different populations of the northern (top) and southern parallels (bottom) as two-dimensional histograms with 15″ bin size overlaid on the DSS2 red greyscale image (see Fig. 1 for the placement of these fields in the 30 Doradus region). Histograms show individual stars selected based on their excess.
In the text
thumbnail Fig. 12

Spatial distribution of different populations in the centre field on top of the F335M mosaic (top left and as background in each panel). The red arrow indicates the location of the central cluster R 136, and the red crosses mark the positions of previously well-studied IR-bright knots (e.g. Walborn et al. 1999). Top right: 12CO moment zero maps from ALMA (data from Wong et al. 2022) as green to yellow contours. The coloured dots show individual stars selected in different populations. The image has a field of view of 7.4′ × 4.4′ (110 × 66 pc).
In the text
thumbnail Fig. 13

Fraction of sources in different selections that overlap with the ALMA CO(2–1) map from Wong et al. (2022). The cross shows the fraction of overlap without any intensity constraint, the plusses show the overlap with regions in which the CO intensity exceeds 1.0 K km s−1. The black symbols refer to the full JWST catalogue, and the grey symbols and bands show the expected values for a uniform distribution of stars across the field of view with both JWST and ALMA coverage. These values are higher than for the old stars (>20 Myr) because only a few are detected in the highly extinguished regions.
In the text
thumbnail Fig. 14

Zoom-ins around the three IR-bright knots 1, 2, and 3. Top row: colour images of 5″ × 5″regions. Red: F444W, green: F200W, blue: F090W. Bottom row: greyscale F200W images with selected sources as coloured dots. The larger purple dots refer to IR-excess sources, and the smaller maroon, red, orange, and yellow dots are stars selected with isochrones with ages of 0.5, 1, 2, and 4 Myr, respectively.
In the text
thumbnail Fig. 15

Colour-magnitude diagram of R 136 (orange dots) on top of the full region (black). The orange dots refer to all stars within 10″ from the centre of R 136. The coloured lines show theoretical MIST isochrones.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.