HUBBLE TARANTULA TREASURY PROJECT: UNRAVELING TARANTULA'S WEB. I. OBSERVATIONAL OVERVIEW AND FIRST RESULTS^*

Starbursts are short-lived periods of intense, massive star formation (Searle et al. 1973) that occur in compact regions (10–1000 pc) and dominate the overall luminosity of their host galaxy (Heckman 2005). In their annual review Kennicutt & Evans (2012) remind us that the star formation rate per unit area in a starburst is much larger than in disk-averaged star formation rate surface densities (∼0.1 M_☉ yr⁻¹ kpc⁻²). Evidence for recent starburst activity is found in 6% of the galaxies in the local universe (z ≪ 1; Lee et al. 2007, 2009) and in 15% of those at high redshift (z > 1; O'Connell 2005; Douglas et al. 2010), indicating that starbursts are a relatively common phenomenon. The impact of a major starburst on the evolution of a galaxy is dramatic, because it will shape the galaxy dynamics, stellar evolution and chemical compositions in ways that are largely dictated by the intensity and the duration of the bursting event (i.e., Dekel & Silk 1986; Mac Low & Ferrara 1999; Kennicutt & Evans 2012).

Starbursts are studied on two complementary fronts. High-redshift surveys collect information over a range of long look-back times for a large sample of objects, but suffer from limited or no spatial information and crude age diagnostics. Studies in the local universe, on the other hand, are affected by the small size of the sample, but can probe discrete structures within the star-forming region using more precise population diagnostics (Hunter & Elmegreen 2006; Tolstoy et al. 2009; McQuinn et al. 2010). In both approaches the characterization of starburst properties is still rudimentary and sometimes contradictory.

Studies in individual nearby dwarf galaxies and in surveys at redshift higher than ∼1.0 indicate that a starburst may last for 5–10 Myr, suggesting that the violent feedback from the massive stars formed during the burst suppresses for a while any further star formation (Tosi et al. 1989; Ferguson & Babul 1998; Annibali et al. 2003, 2009; Stinson et al. 2007). Other studies in the local universe, however, advocate for much longer star-forming episodes (>100 Myr; Konstantopoulos et al. 2009; McQuinn et al. 2012), although it is not clear if their approach allows for sufficiently high temporal resolution. Since the spectral energy distributions (SEDs) of galaxies at high redshift can be interpreted only by assuming a duration for the burst, understanding how and for how long star formation can be sustained during a starburst event has profound implications on our picture for understanding galaxies.

The Tarantula Nebula is by far the most luminous (fHα ∼ 1.3 × 10⁻⁸ erg cm⁻² s⁻¹) and massive star-forming region in the Local Group (Kennicutt & Hodge 1986). Covering an area of ∼40, 000 pc², the Tarantula Nebula is the closest extragalactic giant H ii region and is comparable in size to the unresolved luminous H ii complexes observed in distant galaxies (Oey et al. 2003; Hunt & Hirashita 2009).

With more than 800 spectroscopically confirmed OB stars (Evans et al. 2011), some of which are among the most massive candidates (Crowther et al. 2010), 30 Dor is often described as a mini-starburst (Leitherer 1998). The UV flux coming from 30 Dor is such that it has been proposed as a small scale local analog to the Lyman-break galaxies (Meurer et al. 1997; Shapley et al. 2003; Heckman et al. 2004). Furthermore several authors (e.g., O'Connell & Mangano 1978; O'Connell et al. 1995; Brandl et al. 2004) have suggested that the luminous knots observed in bursting galaxies such as M82 are made of (multiple) 30-Dor-like systems.

By virtue of its location in the LMC (∼50 kpc; Panagia et al. 1991; Pietrzyński et al. 2013), with Hubble we can resolve the Tarantula Nebula into single stars down to the sub-solar mass regime (Andersen et al. 2009; De Marchi et al. 2011). The low inclination angle (∼30°; Nikolaev et al. 2004) limits the line of sight confusion, and the foreground reddening is low because of the high Galactic latitude.

The Tarantula Nebula is a very dynamic region both in terms of runaway-star candidates (e.g., Evans et al. 2011) and gas motion (Chu & Kennicutt 1994). Since it is a multi-stage star-forming region (Walborn & Blades 1997), where loose associations and very dense star clusters of different ages coexist in a relatively small volume, the Tarantula Nebula is a noteworthy window in which to test various scenarios of clustered star formation. In particular, observational evidence shows that over the last ≳ 25 Myr the star-formation history of the region (Figure 1) has been complex:

• 1.  
  Walborn & Blades (1997) identified multiple generations of stars ranging from ⩽0.5 Myr to ∼25 Myr in 1'–2' from the core of 30 Dor, R136.
• 2.  
  Grebel & Chu (2000) derived an age of ∼20–25 Myr for Hodge 301, a cluster ∼3' to the northwest of R136.
• 3.  
  NICMOS observations at the "frontier" between Hodge 301 and R136 showed that this region is strewn with several massive O stars embedded in dense knots of dust, suggesting that stellar feedback from the two clusters may be triggering new episodes of star formation (Brandner et al. 2001).
• 4.  
  ∼7' to the west of R136 the super-bubble created by the ∼10 Myr old OB association LH99 is filled by the expanding supernova remnant N157B, which contains the most energetic pulsar known (Chen et al. 2006).
• 5.  
  Deep HST optical and NIR observations suggest that the core of 30 Dor, R136, is an interacting double cluster (Sabbi et al. 2012).

Zoom In Zoom Out Reset image size

As these previous studies suggest, the Tarantula Nebula region is complicated. While there have been attempts to study in detail some of its specific parts, much of the focus with HST to date has been on its central cluster (Hunter et al. 1995, 1996; Walborn et al. 1999, 2002; Selman et al. 1999; Brandner et al. 2001; Andersen et al. 2009; De Marchi et al. 2011; Sabbi et al. 2012; Selman & Melnick 2013). These studies have yielded significant results, nevertheless a comprehensive understanding of the Tarantula Nebula in its entirety can only be reached by detailed studies that include other star clusters and surrounding field stellar populations.

The ongoing HST panchromatic survey of the Tarantula Nebula HTTP (PI: Sabbi, GO-12939) is built on an existing HST program (PI: Lennon, GO-12499) in filter F775W (∼Sloan Digital Sky Survey i-band). HTTP was awarded 60 orbits of HST time. Figure 2 shows the wavelength coverage of the survey, once it is combined with the GO-12499 dataset.

Zoom In Zoom Out Reset image size

The decision to survey the region over such a wide wavelength baseline was driven by several factors. Adding an optical filter to the F775W dataset allows us to derive CMDs down to m_F775W = 27, that for at the distance of the LMC, for a ∼2 Myr old stellar population with Z = 0.008 corresponds to ∼0.5 M_☉. Our photometry goes 4–5 mag below the oldest turn-off and therefore allows us to reconstruct the star formation history of the entire region and infer the changes in star formation rate with time. We preferred the Advanced Camera for Surveys (ACS) F555W to the wider F606W filter, because the latter, being quite red, would have provided a limited color baseline relative to the F775W filter. Furthermore the Hα line falls in the F606W bandpass and in H ii regions this can affect the accuracy of the photometry.

Supplementing our dataset with two NIR filters allows us to extend our study to the dustiest regions, in order to better constrain the properties of cool stars such as red supergiant, red giant (RGB), and asymptotic giant branch stars. Even more important, only NIR observations can identify regions of ongoing star formation that are still embedded in their dusty cocoons (see Figure 3).

Zoom In Zoom Out Reset image size

Previous NICMOS observations have shown that star formation took place in the recent past at the gaseous interface between R136 and Hodge 301 (Brandner et al. 2001). By extending this study to the entire Tarantula Nebula we can identify pockets of ongoing star formation far from the main clusters. Thus we will test the predicted constructive properties of stellar feedback (e.g., Elmegreen & Lada 1977; Vanhala & Cameron 1998), and evaluate the efficiency of triggered star formation as inferred from modeling (e.g., Dale et al. 2012). We have selected the F110W and F160W filters since they provide the highest throughput with the IR channel.

Observations in the NUV have been designed to infer the temperature of the hotter stars and, once combined with optical measurements, constrain the amount of extinction (Zaritsky 1999; Romaniello et al. 2002). We have chosen the Wide Field Camera 3 (WFC3) filter F336W because of its high throughput and because it is specifically designed to constrain the Balmer break. To better determine the temperature of the hotter stars and to provide the largest baseline possible to characterize the extinction law we added the even bluer filter F275W.

Finally we are observing the Tarantula Nebula with the narrow-band ACS filter F658N (corresponding to Hα). Young (<5 Myr) low-mass pre-main sequence (PMS) stars can be easily identified in the optical CMD as a population of faint red objects, and can be used to trace the extent of young star-forming regions (Sabbi et al. 2007, 2012; Cignoni et al. 2010; Gouliermis et al. 2012). With the aging of a stellar population, however, the colors of the PMS stars become too blue to be distinguished from the low-mass MS stars of the field, an effect that becomes even more severe when extinction is highly variable, as is the case for the Tarantula Nebula. In previous studies (Panagia et al. 2000; De Marchi et al. 2010, 2011) we have shown that a combination of broad- and narrow-band filters can be used to identify low-mass PMS stars with active mass accretion, and that this approach is very efficient in picking up older (>10 Myr) and more diffuse episodes of star formation. In addition, the same strategy will allow us to identify main-sequence and evolved stars with Hα excess, such as Be and B[e] stars (Grebel 1997).

In summary the high sensitivity, spatial resolution, and broad wavelength coverage of HTTP will allow us to reconstruct the Tarantula's star-formation history in space and time on a parsec scale. In particular we will be able to (1) characterize when and where star formation occurred, (2) establish the length and the strength of the star-forming episodes and their spatial scale, (3) depict the life cycle of star clusters, (4) portray the role of stellar feedback in shaping the evolution of a starburst, and (5) probe the universality of the stellar initial mass function. In addition, a Bayesian fitting of the stellar SEDs will allow us to accurately determine R_V, and therefore supply information about the size of the dust grains and their spatial distribution (K. Gordon et al. 2013, in preparation).

The HST program GO-12499 (PI: Lennon) represents the first epoch of a proper-motion survey of massive stars in the 30 Dor region. The survey covers a projected area of ∼14' × 12', corresponding to ∼210 × 180 pc at the distance of the LMC.

A total of 15 pointings were used to map the whole region. All the observations were taken between 2011 October 3 and 29, using WFC3 as primary instrument and ACS in parallel. The orientation angle of the mosaic was chosen to include the very massive runaway VFTS#16 (Evans et al. 2010).

Since the goal of the program is to determine the positions of the stars with exquisite accuracy, all the observations were acquired in a single filter. Filter F775W was chosen in an effort to minimize the effect of the strong and variable nebular emission that pervades the entire region.

Each pointing was observed for one orbit of HST time, and each orbit consisted of one short and four deeper exposures with each camera. The exposure times were dictated by the buffer dumping of WFC3 and by the length of the orbit. To give the survey as much spatial uniformity as possible, observations were stepped in such a way that no star would fall in the inter-chip gap in more than one of the WFC3 exposures.

The full dataset consists of 15 exposures of 35 s, 14 of 500 s, and 45 of 700 s with UVIS and 15 of 32 s, 14 of 377 s, and 45 of 640 s with ACS for a total of 148 exposures. Once stacked together the deep exposures result in a mosaic of ∼28, 000 × 22, 000 pixels, with a 40.00 mas pixel⁻¹ scale. Figure 1 shows the mosaic, with colors added based on ground based data. Figure 4 shows the depth map of the deep stacked image.

Zoom In Zoom Out Reset image size

The CCD detectors on board HST are subject to a flux of energetic particles that progressively damage the silicon lattice of the detectors and create charge traps that redistribute the flux from one pixel to the other during the detector readout process. As a result of this cumulative radiation damage, the charge transfer efficiency (CTE) of the CCDs on board HST is progressively degrading. We used the algorithm described by Anderson & Bedin (2010) to correct the effects of the degrading CTE directly on the images. This algorithm is automatically applied by the ACS calibration pipeline. The routine we used to correct WFC3 data for CTE can be downloaded from the WFC3 Web site.¹⁹

As was done for program GO-12499, we are observing the field with WFC3 as primary instrument and ACS in parallel. Visits are organized in four blocks of visits. Figure 5 shows the 36 NIR and 24 UVIS WFC3 pointings superimposed on the DSS image of 30 Dor, together with the ACS parallels. Each pointing corresponds to a single orbit of HST time and consists of a four-step dither-pattern.

Zoom In Zoom Out Reset image size

The first and third blocks are observed with WFC3/IR and ACS. During each orbit we acquire two exposures in the WFC3 F110W and two in the F160W filters, while ACS is used to collect four observations in the F658N filter (in Figure 5, upper panels, WFC3 is in green and ACS is in magenta).

In the second and fourth blocks, each orbit begins with the acquisition of a short exposure in the F336W filter with WFC3 and a short ACS image in the F555W filter. We then collect two images in the F336W and two in the F275W filters. At the same time, we acquire four exposures in the F555W filter with ACS (in Figure 5, lower panels, WFC3 is in blue and ACS is in magenta).

Half of the IR observations were acquired between 2012 December 12 and 19 and are described in the next section. The remaining half is scheduled for the summer of 2013. The NUV observations will be collected in 2013 April and September. This observing strategy is designed to exploit the natural rotation of the telescope during the year.

The filters, the number of exposures, and the exposure times used in HTTP are summarized in Table 1. Updates on the status of the survey can be found at http://30dor.stsci.edu.

Observations with very low background (<10 e⁻ for WFC3 and <50 e⁻ for ACS (for a discussion of the background levels in UVIS observations, see Baggett & Anderson 2012), such as those acquired in the NUV and/or through narrow-band filters) suffer large losses for very faint sources, to the point that faint sources can completely vanish during the readout transfers. Since Cycle 20, it has been possible to mitigate the effects of the degrading CTE by adding an extra amount of photons through an internal illumination. This process, called "post-flash," has the effect of artificially increasing the background level (MacKenty & Smith 2012). All the NUV observations of HTTP are acquired using "post-flash." We decided to not use post-flash for the Hα exposures, since we estimated that the background would be sufficiently high. Even with the background increased, all sources still suffer some CTE loss, and it is necessary to make a photometric correction for the detected sources. For this reason we will correct all the data acquired in the NUV and visual regimes using the Anderson & Bedin algorithm.

Observations in the F775W filters cover a projected area of ∼14' × 12'. The first step in reducing the data was to create a distortion-free reference frame and relate the astrometry and photometry of each exposure to this frame. To achieve this, in each F775W deep exposure we measured positions and fluxes for all the sources that had more than 100 counts in each exposure and no brighter neighbors within a 5 pixels radius using the program img2xym_WFC.09x10 (Anderson & King 2006).²⁰ This program uses libraries of empirical point-spread functions (PSFs) to account for the spatial variations caused by the optics of the telescope and the variable charge diffusion in the CCD and creates a geometric-distortion corrected list of sources.

We used the Two Micron All Sky Survey (2MASS) catalog (Skrutskie et al. 2006) as an initial astrometric reference frame. We then found common stars between each image list and 2MASS to define the general six-parameter linear transformation between the distortion-corrected frame of each exposure into the 2MASS-defined frame. We used these transformations to collate the individual star lists in the reference frame and determine the average position for ∼110, 000 bright, isolated, unsaturated stars that in the deep exposures have three or more coincident detections.

For the ACS data we initially used the geometric distortion and PSF library of Anderson & King (2006). Both the PSF and distortion solution were determined using data acquired before the last Hubble servicing mission (SM4). A careful inspection of the position residuals showed that the PSF has developed a 2% asymmetry and the distortion solution has changed by ∼4% of a pixel since SM4. Although this is a very small change, it is sufficient to introduce noticeable errors in very wide mosaics, therefore we used all the bright and isolated stars found in the ACS field of view to derive a new PSF and update the geometric distortion. Since all the observations in this program were acquired with very similar rotation angle, we were not able to independently solve for the ACS skew term, however, we used the common stars between ACS and WFC3 to calibrate it relative to WFC3.

Finally we improved the internal quality of the reference frame by iteration. The rms residuals of the average positions were less than 0.01 pixel in each coordinate. In this process we have mapped the ACS pixel scale to the size of the WFC3 UVIS channel, so that in all the catalogs a pixel corresponds to 0040. The final reference list was then used to compute the final astrometric transformation of each exposure (for both ACS and WFC3) into the reference frame.

The stars used in Section 4.1 to align the single exposures to the reference frame were found in a single pass through each image. This approach will find almost all the bright stars in a field, but will miss many of the obvious faint sources, especially those in the wings of bright objects.

To identify the fainter sources and push the photometry as deep as possible, we used an evolution of the program described in Anderson et al. (2008). The new program will be fully described in a forthcoming paper (J. Anderson et al., in preparation).

The previous version of the code was designed to work only with ACS and could handle only two filters. This new routine can handle multiple instruments and several filters at the same time. The routine first finds the bright stars, then subtracts them and searches for fainters stars in the subtracted images. It performs several such finding iterations (Figure 6), using detection thresholds that depend on the number of exposures available at each point in the field. This approach is particularly useful when, as in the case of our survey, the coverage is not uniform (in some regions we have only one or two exposures, while in others we have up to 14 overlapping exposures, Figure 4). An inspection by eye shows that the code efficiently avoids spurious detections caused by the diffraction spikes of saturated stars and sharp gaseous filaments (Figure 6).

Zoom In Zoom Out Reset image size

We were able to identify more than 1.1 × 10⁶ sources in the F775W observations. The photometry was calibrated to the VEGAmag photometric system using the zero points listed on STScI Web site.²¹ These zero points are derived for images combined using "astrodrizzle," thus they cannot be directly applied to our flt-based photometric catalogs. To calibrate the flt-based catalogs into the VEGAmag system we selected several isolated bright stars in the drizzled images and measured their magnitudes using 05 aperture photometry for ACS and 04 for WFC3 data and applied the appropriate zero point. For these stars we then determined the difference between the flt-based magnitude and the magnitude in the VEGAmag system and used the average differences as the final zero points.

The F775W photometric catalog was used to measure the sources in the IR images. Of the initial 1.1 × 10⁶ sources, ∼480, 000 are in the region observed with the first block of IR data.

In addition to finding and measuring stars, the reduction routine also produces several diagnostics to help identify the stars that are well-measured. One such parameter is Q, which reports the linear-correlation coefficient between the PSF and the stellar image: Q = 1.0 is a perfect fit. This parameter is particularly useful in identifying unresolved objects such as background galaxies and blended sources. We also used the rms of the individual photometric measurements about the mean to help identify well measured stars.

In both ACS and WFC3 images acquired with the F775W filters, stars brighter than m_F775W = 15 are saturated even in the shorter exposures. As a result the Q parameter for these stars differs significantly from 1. We used aperture photometry to recover the flux of these sources, and we did not apply any selection criteria. For ACS we then selected the sources in the magnitude range between $15\le m_{\rm F775W_{ACS}} <21$ with Q > 0.95, while for sources fainter than $m_{\rm F775W_{ACS}}\ge 21$ we selected only those objects with Q > 0.92 or photometric error <0.5. Similarly, for the F775W WFC3 filter in the magnitude range between $15\le m_{\rm F775W_{WFC3}} <22$ we selected only the objects with Q > 0.95 and for fainter sources the photometric error had to be better then 0.5.

Because of the coarser resolution of the IR channel, crowding is more severe in the filters F110W and F160W and the parameter Q does not provide an estimate of the quality of the photometry that is as good as in the case of the optical wavelengths. Moreover, for most of the field in the NIR we have only two measurements per filter, therefore to estimate the photometric error we used the signal to noise instead of the rms of the individual photometric measurements. As a result for both the F110W and F160W filters we selected only those objects with Q > 0.85 and photometric error smaller than 0.1. After these selections the combined ACS+WFC3 catalog in the F775W filter contains more than 366,000 stars, of which ∼150, 000 that are detected in both the F110W and F160W filters. Figure 7 shows Q as a function of magnitude for the two F775W filters (ACS and WFC3 are showed separately) and for the NIR data. Sources that meet the selection criteria are shown in black, while the remaining sources are in gray.

Zoom In Zoom Out Reset image size

To better understand the spatial distribution of the identified stellar populations, and at the same time, highlight how the reddening changes as a function of position, we have divided our catalog into 18 regions of 3000 × 3000 pixel each, corresponding to ∼29 pc × 29 pc. Figure 9 shows the location of these regions on the F160W mosaic (upper panel) and the corresponding 18 m_F110W versus m_F110W − m_F160W CMDs (lower panel).

The inspection of these CMDs shows that the UMS is always well populated, with the exception of regions A, B, D, G, Q, and R. These six regions are likely sampling the field of the LMC.

The majority of R136 is in region I, and the NE clump identified by Sabbi et al. (2012) is in region K. The CMDs for both these regions show a large number of red and faint PMS stars, as well as very well populated UMSs. Both R136 and the NE clump are very young systems (⩽5 Myr), and the PMS stars are well separated from the LMS stars. PMS stars are likely also present in regions L and N.

Hodge 301 is in region O. The UMS of the corresponding CMD is well populated in agreement with the young age of the cluster. At these ages PMS stars are too blue to be distinguished from the LMS stars using broad-band photometry only.

The RC and the RGB vary from relatively compact and well defined (as in regions B, E, M, O, Q, and R) to extremely broad. For example in region J, the RC covers more than 4 mag in the F110W filter, a typical signature of differential reddening. In their analysis of the reddening in the LMC, based on the RC, Haschke et al. (2011) find that not only 30 Dor exhibits the highest reddening of any region in this galaxy, but also the highest differential reddening.

It is interesting to note that young and old stars are affected by different amounts of extinction, with the reddening being on average less severe for the young stars, than the old ones. This effect is particularly evident for example in region H, where the UMS is narrow and well defined, but the magnitude of the RC ranges between 18 ≲ m_F110W ≲ 20. This is apparently at variance with conclusions by Zaritsky (1999), who found that on average in the LMC old stars are less obscured by dust than hot young stars.

In deriving the reddening map, Zaritsky (1999) used only those stars that were detected at the same time in the four bands U, B, V, and I. This selection criterion penalizes the most extinguished evolved stars, therefore the average amount of reddening derived for the RGB stars may have been underestimated. As shown in Haschke et al. (2011) in these highly extinguished regions the distribution of the reddening values is not Gaussian, but shows an extended tail toward the higher values.

It has to be noted that most of the UMS stars found in the Tarantula Nebula region are associated with the starburst that originated R136 and can be all considered to be at the same distance from us. The RGB stars, on the other hand, are evenly distributed in the entire body of the LMC and are more affected by the three-dimensional structure of the galaxy. If for example 30 Dor is on the nearer side of the LMC disk, the majority of the RGB stars in our field would be beyond 30 Dor and will be extinguished by the dust associated to the Tarantula Nebula. We will be able to better characterize this issue when we observe the region in bluer filters (F275W, F336W, and F555W).

The CMDs in Figure 9 suggest that stellar populations of different ages are not uniformly distributed. To better highlight this effect, in Figure 10 we have compared the spatial distribution of stars selected from different evolutionary sequences.

Zoom In Zoom Out Reset image size

In panel (A) we have selected UMS stars brighter than m_F110W < 16.6, and bluer than m_110W − m_160W < 0.3. These stars are younger than ∼20 Myr and are mostly concentrated around R136 and the NE clump. The IR observations collected so far cover 9 (namely, 537, 470, 398, 385, 382, 392, 706, 682, and 581) of the 15 O stars identified by Bressert et al. (2012) as candidate high-mass stars formed in isolation, and another 5 targets will be observed in 2013 July. As already reported by these authors, there is no evidence for an underlying cluster around these sources. Further constraints on the true nature of these stars will come when proper motion measurements will be available.

Panel (B) shows the distribution of stars in the magnitude range between 17.9 < m_F110W < 19.9, and bluer than m_F110W − m_F160W < 0.4. At the distance of the LMC this correspond to a mass range between ∼3 and ∼7 M_☉. After ∼2 Myr the most massive PMS stars in this mass range are already merging into the UMS, while a ∼3 M_☉ star will still be on the MS after ∼150 Myr (Dell'Omodarme et al. 2012). Therefore panel (B) shows the distribution of stars in the age range between ∼2 and 150 Myr. R136 and the NE clump are still very well defined. The cluster to the northwest of the image is Hodge 301. A careful inspection of panel (A) shows that Hodge 301 is also visible in that map.

Sources in panel (C) have colors between 0.7 < m_F110W − m_F160W < 1.3 and magnitudes between 20.5 < m_F110W < 22.5. This is the locus of ≲ 1 M_☉ PMS stars. These sources are likely younger than ∼5 Myr. Hodge 301 is not visible in this region, while a "chain" of young clustered objects connects the NE clump to the northern corner of the map. Some of these clusters coincide with the embedded massive O stars found in the dense knots of dust at the interface between R136 and Hodge 301, already identified by Brandner et al. (2001). Others are aligned with one side of the shell of ionized gas that surrounds a soft and bright X-ray-emitting bubble (Meaburn 1984; Wang & Helfand 1991; Townsley et al. 2006). Some of these clusters likely host embedded OB stars (Walborn et al. 2013). Although it is possible that the feedback from the two massive clusters NGC 2070 and Hodge 301 may have triggered the star formation in the nearer clumps, the majority of the clusterings are at a projected distance from R136 between 30 and 60 pc, too far away to be affected by R136 in the less then ∼2 Myr. Alternatively, these systems may result from the collapse of residual pockets of gas formed during the fragmentation of the giant molecular cloud that formed R136, as predicted by the models that invokes fractal initial conditions and assembling of sub-clusters to form massive young clusters (Aarseth & Hills 1972; McMillan et al. 2007; Allison et al. 2009; Moeckel & Bonnell 2009; Bonnell et al. 2011; Smith et al. 2011; Fuji et al. 2012; Fujii & Portegies Zwart 2013).

LMS stars (panel (D)) have been selected in the color range 0.26 < m_F110W − m_F160W < 0.65 and magnitudes between 21.5 < m_F110W < 22.5. These sources have ages between a few tens of megayears and several gigayears. The spatial distribution of these sources is quite uniform in agreement with our assumption that we are looking at the stars in the field of the LMC. The only visible cluster in panel (D) is Hodge 301, in agreement with the fact that old PMS stars cannot be separated from LMS stars using broad-band photometry only. R136 and the NE clump are not visible in this map, as expected from their very young (<5 Myr) age.

The number of LMS stars is significantly lower in the regions that in panel (C) correspond to sites of very recent star formation. These regions are characterized by a higher stellar density and nebulosity, making it harder to detect faint stars, and thus the low number of LMS may be the result of a lower level of completeness. We note that many of these regions are obscured by dust. If the Tarantula Nebula is on the nearer side of the LMC disk, the dust associated to the star-forming region would attenuate the luminosity of the LMS, pushing them outside the range of magnitudes and colors that we are considering. Artificial star tests and stellar SED fitting will allow us to better discriminate between these two hypotheses. The spatial distribution of RGB stars (panel (E), 0.6 < m_F110W − m_F160W < 1.1, 15.3 < m_F110W < 18.5) is noisy because of low number statistics, but on average appears quite uniform.