The VMC Survey. XXII. Hierarchical Star Formation in the 30 Doradus-N158–N159–N160 Star-forming Complex

It has been suggested that star formation is hierarchical, as revealed by stellar structures spanning a wide range of scales, from large stellar complexes and aggregates to small associations and clusters (Efremov 1995; Elmegreen et al. 2000). Hierarchical stellar structures have been investigated for a few famous star-forming regions, such as Taurus, Orion, Ophiuchus, etc. in the Milky Way (Gomez et al. 1993; Larson 1995; Simon 1997) and NGC 345 in the Small Magellanic Cloud (SMC; Schmeja et al. 2009; Gouliermis et al. 2014). On larger scales, a number of galaxies display galaxy-wide hierarchies formed by young stars (e.g., Elmegreen & Elmegreen 2001; Elmegreen et al. 2006, 2014; Bastian et al. 2007, 2009; Gieles et al. 2008; Gouliermis et al. 2010, 2015). These studies show that hierarchical stellar structures, significant for a high degree of substructure and an absence of characteristic scales, exhibit self-similar and fractal properties. These features are also found for the structures formed by the interstellar medium (ISM), such as clumps and filaments on various scales (e.g., Elmegreen & Falgarone 1996; Roman-Duval et al. 2010). This similarity suggests that the hierarchical structures of the young stars may originate from those of their parental ISM, which are in turn related to turbulence (Larson 1981; Elmegreen & Scalo 2004; Federrath et al. 2009; Girichidis et al. 2012), agglomeration with fragmentation (Carlberg & Pudritz 1990; McLaughlin & Pudritz 1996), and self-gravity (de Vega et al. 1996; see also the review of Elmegreen et al. 2000). After birth, the hierarchical stellar structures evolve toward uniform distributions, with the substructures becoming eliminated in roughly one crossing time (Gieles et al. 2008; Bastian et al. 2009; Gouliermis et al. 2015).

30 Doradus (30 Dor hereafter, also known as the Tarantula Nebula) in the Large Magellanic Cloud (LMC), with R.A. (J2000) = 05^h38^m38^s, decl. (J2000) = −69^o05'42'', is among the most famous astronomical objects. It is by far the most luminous and massive star-forming region in the Local Group (Kennicutt & Hodge 1986) and hosts at its center the star cluster NGC 2070, which, in turn, is a collective of dense subclusters (Walborn & Blades 1997; Sabbi et al. 2012). On the other hand, three active star-forming regions, N158, N159, and N160, are found to the south of 30 Dor (Nakajima et al. 2005; Galametz et al. 2013). They are located in the northern part of the LMC's remarkable "molecular ridge," which stretches ∼2 kpc in a nearly north–south straight line and contains almost one third of the LMC's molecular gas traced by CO (Cohen et al. 1988; Kutner et al. 1997; Johansson et al. 1998; Fukui et al. 2008; Ott et al. 2008). Detailed studies of the hierarchical stellar structures in these star-forming regions have not been reported. Although Bastian et al. (2009) have investigated this issue for the entire LMC, their stellar sample did not reach the necessary resolution to reveal the substructures of these regions. This paper thus aims to carry out such a study for the 30 Dor-N158–N159–N160 star-forming complex based on young upper main-sequence (UMS) stars, which have been carefully selected from the VISTA Survey of the Magellanic Clouds (VMC; Cioni et al. 2011). We use dendrograms to illustrate the "parent–child" relations of star groups and subgroups identified on various density levels; and, additionally, the group size distribution is explored. This will demonstrate the hierarchical structures formed by the young UMS stars, which may indicate a scenario of hierarchical star formation in this complex.

This paper is organized as follows. Section 2 describes the VMC survey and the data used in this work. Selection of UMS stars is outlined in Section 3, while in Section 4 we show their spatial distributions. We identify groups of UMS stars and reveal their "parent–child" relations using dendrograms in Section 5, and the group sizes are investigated in Section 6. We finally close this paper with a summary and conclusions.

Data used in this work come from the VMC survey (Cioni et al. 2011), which is carried out with the 4.1 m Visible and Infrared Survey Telescope for Astronomy (VISTA; Sutherland et al. 2015). The VMC survey is a multi-epoch, uniform, and homogeneous photometric survey of the Magellanic system performed in the near-infrared Y, J, and ${K}_{{\rm{s}}}$ filters (centered at 1.02 μm, 1.25 μm, and 2.15 μm, respectively). On completion, the survey is expected to cover ∼170 deg² of the Magellanic system. A sequence of six offset positions, each one called a pawprint, is used to fill in the gaps between the 16 detectors of the VISTA infrared camera (VIRCAM; Dalton et al. 2006); the combined image, or a VISTA tile, covers an area of ∼1.5 deg². Each tile in the VMC survey is designed to be observed at 3 epochs in the Y and J bands, and 12 epochs in the ${K}_{{\rm{s}}}$ band. The total exposure times are 2400 s in the Y and J bands, and 9000 s in the ${K}_{{\rm{s}}}$ band for most regions of each tile, except for the tile edges and some areas of extra overlap among the detectors. The detectors usually have saturation limits of Y = 12.9 mag, J = 12.7 mag, and ${K}_{{\rm{s}}}=11.4$ mag, and the stacked images from all epochs can provide sources with typical 5σ limiting magnitudes of Y = 21.9 mag, J = 22.0 mag, and ${K}_{{\rm{s}}}=21.5$ mag; but the saturation limits and photometric depths also depend on sky conditions and crowding (Cioni et al. 2011; Tatton et al. 2013). This depth reaches the oldest main-sequence (MS) turnoffs in both the LMC and SMC (Kerber et al. 2009).

The 30 Dor-N158–N159–N160 star-forming complex is analyzed based on tile LMC 6_6 (see Cioni et al. 2011 for tile definitions). We retrieved the data of tile LMC 6_6 as part of VMC Data Release 3 from the VISTA Science Archive (VSA). The details of the VSA and the VISTA data flow pipeline can be found in Cross et al. (2012) and Irwin et al. (2004), respectively. We have used point-spread function (PSF) photometry from Rubele et al. (2012) instead of aperture photometry to reduce the influence of source crowding (Tatton et al. 2013). Rubele et al. (2012) generated deep, PSF-homogenized tiles by stacking images from different epochs, each one formed by the combination of six pawprints. They performed PSF photometry on the deep tiles and estimated the photometric errors and local completeness using artificial star tests. Specifically, the local completeness is estimated in a ring of radius 0025 around each source and in bins of ±0.05 mag separately in the Y, J, and ${K}_{{\rm{s}}}$ bands. We use their local completeness estimates to assign weights to stars in our sample (see Sections 3 and 4). A detailed description of the PSF homogeneity and photometry can also be found in Rubele et al. (2015).

The top half of VISTA's detector 16 has varying quantum efficiency on short timescales, making accurate flat fielding impossible. This leads to worse signal-to-noise ratios and unreliable data in regions covered by that detector (Rubele et al. 2012; Tatton et al. 2013). For tile LMC 6_6, the region affected is the southwest corner with R.A. (J2000) < 05^h33^m55^s and decl. (J2000) < −69^o43'48''. This region is excluded from our analysis.

The left-hand panel of Figure 1 shows the ($J-{K}_{{\rm{s}}}$, ${K}_{{\rm{s}}}$) CMD of sources in tile LMC 6_6 without correcting for extinction. The MS, red-giant branch (RGB), and red clump (RC) populations are all resolved. We use the ($J-{K}_{{\rm{s}}}$) color instead of ($Y-{K}_{{\rm{s}}}$) because the former suffers less from interstellar extinction. Still, effects of interstellar extinction varying across the tile are present, evidenced by the width of the MS and the RGB, and the elongated tail associated with the RC toward fainter magnitudes and redder colors.

Zoom In Zoom Out Reset image size

We estimate the extinction by taking advantage of the extinction values of ∼1.5 × 10⁵ RC stars in this tile provided by Tatton et al. (2013). We bin both all stars and their RC stars into spatial grids with a grid size of 1' × 1'; the median extinction value of RC stars is taken as the extinction estimate for all stars in each bin. We use the extinction coefficients A_J = 0.283 $\times \,{A}_{V}$ and ${A}_{{K}_{{\rm{s}}}}=0.114\times {A}_{V}$, computed from the Cardelli et al. (1989) extinction curve with ${R}_{{\rm{V}}}$ = 3.1 (Girardi et al. 2008). De Marchi et al. (2016) reported a new extinction law for 30 Dor, characterized by R_V ∼ 4.5. However, using this law would not lead to significant changes in the corrections for the near-infrared wavelengths considered in this paper (Gordon et al. 2003). Contamination by RGB stars is not important, since their intrinsic color exhibits only a small difference of ∼0.05 ± 0.02 mag from that of the RC stars (Tatton et al. 2013), thus introducing very minor, if any, errors in correcting the CMD. The middle panel of Figure 1 shows the CMD after correcting for extinction. The MS and RGB become tighter and the elongated tail associated with the RC appears less prominent.

The extinction in the LMC is subject to a population dependence, which is, however, difficult to model (Zaritsky et al. 2004; Cignoni et al. 2015). Compared with intermediate-age RC stars, the young stars form in gas-rich regions with high extinction, but their ionizing fluxes and stellar winds tend to reduce extinction by evacuating the surrounding material. The problem is further compounded by their different spatial distributions, with respect to the dust. Along the line of sight, dust may lie in the foreground and/or background of a star, and only the foreground dust contributes to its extinction. Perpendicular to the line of sight, dust in star-forming regions also exhibits a high degree of spatial variability (e.g., Lombardi et al. 2010, 2011; Cambrésy et al. 2013). This is also demonstrated by the Herschel dust emission map across the SMC star-forming region NGC 345 (also known as N66), which is highly variable on very small scales and rich in structures of knots, arcs and filaments, etc. (Hony et al. 2015). It is very challenging and beyond the scope of this work to take into account all of these considerations by obtaining a population-dependent, three-dimensional extinction map with high angular resolution. Still, our treatment of extinction is reasonable as the features become more well-defined in the extinction-corrected CMD.

In the right-hand panel of Figure 1, we overplot PARSEC isochrones (version 1.2S; Bressan et al. 2012) of metallicity Z = 0.008 and ages log(τ/yr) = 6.6, 7.0, 7.5, 8.0, and 9.0, shifted by an LMC distance modulus of (m − M)₀ = 18.49 mag (Pietrzyński et al. 2013; de Grijs et al. 2014) and assuming zero extinction. We have subtracted offsets of 0.026 mag in J and 0.003 mag in ${K}_{{\rm{s}}}$, given by Rubele et al. (2015), from the isochrones to correct for the small differences between the VISTA system and the model Vega system. The differences are caused by the unaccounted for nonlinearities in the color–color relation used to calibrate VISTA's photometric zeropoint with 2MASS stars (for details, see Rubele et al. 2012, 2015). It is clear that the upper part of the MS corresponds to a very young population, since it coincides well with isochrones of log(τ/yr) = 6.6 and 7.0. We thus select a sample of UMS stars using color and magnitude cuts, −0.45 ≤$\,J-{K}_{{\rm{s}}}\leqslant 0.05$ mag and $12.5\leqslant {K}_{{\rm{s}}}\leqslant 17.5$ mag, which are shown as the blue solid-line box in the right-hand panel of Figure 1. This sample contains ∼1.5 × 10⁴ stars.

We confirm that the UMS stars in our sample comprise indeed a young population by comparing their ${K}_{{\rm{s}}}$-band luminosity function with synthetic ones (Figure 2). For the observed UMS sample, we take into account the photometric completeness by assigning to each star a weight w = 1/f, where f is the local completeness; because we are using the J and ${K}_{{\rm{s}}}$ bands, f is estimated from the lowest value of the J- and ${K}_{{\rm{s}}}$-band local completenesses, f_J and ${f}_{{K}_{{\rm{s}}}}$, as estimated by Rubele et al. (2012, 2015; Section 2). On the other hand, we create four synthetic populations of metallicity Z = 0.008 and ages log(τ/yr) = 6.6, 7.0, 7.5, and 8.0. We assume a Chabrier (2001) stellar initial mass function (IMF), a 30% binary fraction, and that the binaries are non-interacting systems with primary/secondary mass ratios evenly distributed from 0.7 to 1.0 (Kerber et al. 2009; Rubele et al. 2012, 2015). In each synthetic population, sufficiently large numbers of stars are created and placed at the LMC's distance modulus of 18.49 mag (Pietrzyński et al. 2013; de Grijs et al. 2014). Synthetic UMS samples are obtained by applying the same selection criterion to the synthetic populations as to the real data; their luminosity functions are then calculated.

Zoom In Zoom Out Reset image size

Figure 2 shows that the luminosity functions of the log(τ/yr) = 6.6 synthetic UMS sample and the observed one are consistent, although the latter seems to contain slightly more stars with ${K}_{{\rm{s}}}$ < 16.0 mag (Figure 2(a)); the log(τ/yr) = 7.0 synthetic sample also has a very similar luminosity function, except for a slight lack of stars brighter than ${K}_{{\rm{s}}}=14.0$ mag (Figure 2(b)). Synthetic samples with log(τ/yr) = 7.5 and 8.0 deviate from the observed sample even more significantly, in that they have many fewer stars brighter than ${K}_{{\rm{s}}}=15.7$ mag or 16.8 mag, respectively (Figures 2(c) and (d)). The comparison leads us to suggest that the observed UMS sample is indeed composed of very young stars with ages of a few million years.

The VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011) has collected fiber spectroscopy of about 800 early-type stars in a 25' diameter field of the 30 Dor region, and spectroscopic classification of 352 O–B0 and 438 B-type stars has been done by Walborn et al. (2014) and Evans et al. (2015), respectively. We compare our UMS sample with their OB-type stars; 279 O-type and 313 B-type stars observed by VFTS are covered by our UMS sample, and their spectral types range from O2 (earliest) to B5 (latest). This comparison again indicates the youth of the UMS sample.

Contamination by foreground Galactic stars and background galaxies is minimal. The former reside at typically ($J-{K}_{{\rm{s}}}$) = 0.7 mag, and background galaxies have colors even redder than this value. Thus, they do not overlap with the UMS stars in the CMD and are not included by the color–magnitude selection criterion.

Schmeja (2011) compared different algorithms in constructing stellar density maps and identifying star clusters as overdensities; he showed that clusters with low overdensities or highly hierarchical structure are only reliably identified by methods with inherent smoothing, such as simple star counts and the $n\mathrm{th}$ nearest-neighbor method (Casertano & Hut 1985; Megeath et al. 2016). Thus, in this section, we apply these two methods to the UMS sample to estimate their surface densities; the surface-density map derived with the $n\mathrm{th}$ nearest-neighbor method is then used to identify stellar groups in Section 5. In the first method, we use a spatial grid of $1^{\prime} \,\times$ 1' cells and determine the total weight of stars in each bin, with each star assigned a weight in the same way as described in Section 3. In the second method, surface densities are calculated around each star with

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{n}=\displaystyle \frac{{w}_{\mathrm{tot}}}{\pi {r}_{n}^{2}},\end{eqnarray} \tag{ 1 }$$

where r_n is the distance to the $n\mathrm{th}$ nearest neighbor, and ${w}_{\mathrm{tot}}$ is the sum of the weights of all n+1 UMS stars at a radius of $r\leqslant {r}_{n}$:

$$\begin{eqnarray}&&{w}_{\mathrm{tot}}=({w}_{0}-1)+\displaystyle \sum _{i=1}^{n-1}\,{w}_{i}+0.5\times ({w}_{n}-1).\end{eqnarray} \tag{ 2 }$$

Casertano & Hut (1985) argued that the $n\mathrm{th}$ nearest-neighbor method provides an unbiased estimate of density if neither the central star nor its $n\mathrm{th}$ nearest neighbor are counted. In other words, we take into account the total weight inside the annulus between the central star and its $n\mathrm{th}$ nearest neighbor. Thus, unity is subtracted from the central star, and its excess weight above unity, i.e., ${w}_{0}-1$, represents the over-density inside the annulus; similarly, 1 is subtracted from its $n\mathrm{th}$ nearest neighbor, and half of its excess weight above 1, i.e., 0.5 × (${w}_{n}-1$), is assumed to be inside the annulus. We use n = 10, which has been shown to be a reasonable choice to balance statistical noise and locality (Gouliermis et al. 2010; Megeath et al. 2016). The separations between stars and their tenth nearest neighbors range from 01 in the densest area, to 36 in low-surface-density regions, with a median value of 09. Surface-density maps derived with the two methods are shown in Figure 3.

Zoom In Zoom Out Reset image size

Figure 3 shows that the surface densities are generally low in the outer part of the tile. The UMS stars in the low-surface-density areas could come from the LMC's general background field, which are not associated with the star-forming complex. As we will show later in this paper, the young stars are hierarchically grouped in a way that small and compact stellar structures reside in larger and looser ones. Thus, there is also a possibility that the UMS stars in the low-surface-density areas belong to a larger stellar structure created in the star-forming process, but our sky coverage is not large enough to reveal it. In this work we do not try to reach scales larger than a single tile, thus we characterize this low surface density only from the LMC's background contribution. This should not make significant differences to our results, considering its very small values compared with the higher surface densities of the star-forming complex. We estimate the background level using four subregions, which are devoid of surface-density enhancements and outlined with white boxes in the left-hand panel of Figure 3. We calculate the average and standard deviation of the surface densities derived with simple star counts of the bins in the subregions, iteratively rejecting 3σ outliers. The average value is log($\overline{{\rm{\Sigma }}}\cdot$pc²) = −2.30 and the standard deviation is log($\sigma {\rm{\Sigma }}\cdot$pc²) = −2.28. The former is used as a measure of the background level, and the latter is an estimate of its fluctuations.

High surface densities are found in the central part of the tile, with roughly $-20^{\prime} \leqslant {\rm{\Delta }}$R.A. ≤ 25', $20^{\prime} \leqslant {\rm{\Delta }}$decl. ≤ 30'. The peak surface density, which occurs in the center of 30 Dor with ΔR.A. = 6' and Δdecl. = 17', reaches log(${\rm{\Sigma }}\cdot$pc²) ∼ −0.20, derived from simple star counts, or log(${\rm{\Sigma }}\cdot$pc²) ∼ 0.25, obtained with the $n\mathrm{th}$ nearest-neighbor method, where n = 10. The discrepancies between results obtained with the two methods arise from their different smoothing scales. The star-count method has a fixed smoothing scale of 1'; the $n\mathrm{th}$ nearest-neighbor method, in contrast, has variable smoothing scales, which are the separations between stars and their $n\mathrm{th}$ nearest neighbors. As previously mentioned, the separations range from 01 to 36 for n = 10, depending on the surface densities.

Also note that the majority of UMS stars are distributed in a number of highly fractured groups. The largest ones, with scales of several to ∼20', are consistent with the OB associations cataloged by Lucke & Hodge (1970), whose designations have been labeled in the figure. Here we recall that it can be difficult to define clusters, associations, or larger structures, and equally difficult to make clear distinctions between them (Bressert et al. 2010; Krumholz 2014). Gieles & Portegies Zwart (2011) suggested that bound and unbound stellar systems can be distinguished by comparing the ages and crossing times; yet observationally, it is not trivial to obtain these parameters. In order to avoid ambiguity, we will refer to any stellar structure with a surface-density enhancement as a group, regardless of its size or dynamical state; the critical surface density to define an "enhancement" is not fixed but can be varied arbitrarily, and by doing this we will reveal the hierarchy of stellar groups on a range of scales in Sections 5 and 6. But here we temporarily set the critical surface density at 3σ above the average of the background level, so that we can investigate the stellar distributions on large scales. This corresponds to log(${{\rm{\Sigma }}}_{\mathrm{cr}}\cdot$pc²) ∼ −1.7, which is shown as white or black contours in both panels of Figure 3.

Most prominently in Figure 3, the 30 Dor Nebula contains three populous groups, LH 100, LH 99, and LH 90, from its center to the west. In the vicinity of 30 Dor, there are three large, diffuse stellar groups, LH 89, LH 111, and LH 96, to its northwest, east, and south, respectively. In the southeastern part of the tile, there are another few groups with significant stellar overdensities, i.e., LH 101, LH 104, LH 103, and LH 105, which correspond to the N158, N160, and N159 star-forming regions (Bolatto et al. 2000; Nakajima et al. 2005; Galametz et al. 2013). A number of smaller groups can also be seen across the tile, e.g., LH 74 in the west, and even more that have not been cataloged by Lucke & Hodge (1970).

Many of the groups described above are, in turn, highly subgrouped. For example, LH 96 is far from centrally concentrated and contains multiple surface-density peaks. This has been noted by Lucke & Hodge (1970), who identify four smaller associations inside, i.e., LH 93, LH 94, LH 97, and LH 98. Similarly, the star-forming region, N158, contains two surface-density peaks, LH 101 and LH 104, as cataloged by Lucke & Hodge (1970) and also revealed here. Although not obvious in Figure 3, LH 100 contains two star clusters, NGC 2070 at the center of 30 Dor and Hodge 301 slightly to the northwest (Cignoni et al. 2015); the former is, in turn, a collective of smaller star clusters (Walborn & Blades 1997; Sabbi et al. 2012). A statistical analysis of the hierarchical structures is given below.

In this section, we investigate the hierarchy of UMS star groups through the so-called dendrograms (Houlahan & Scalo 1992; Rosolowsky et al. 2008; Gouliermis et al. 2010). Dendrograms are structure trees showing the "parent–child" relations of groups with surface-density enhancements over different critical values. We vary the critical surface density from log(${{\rm{\Sigma }}}_{\mathrm{cr}}\cdot$pc²) = −1.7 to −0.1; the lower limit corresponds to 3σ above average of the background level, and the upper limit is chosen close to the peak surface density (see Section 4). Between them we define levels of critical surface densities with equal logarithmic steps of 0.2 dex. This step is arbitrarily chosen in order to avoid too many levels in displaying the dendrograms, considering that the surface densities span two orders of magnitude; reducing this step and increasing the number of levels would not change the conclusion reached in this section. On each level we use a friends-of-friends group identification algorithm (Megeath et al. 2016). To identify stellar groups, we use the surface densities estimated with the $n\mathrm{th}$ nearest-neighbor method (n = 10; Section 4). For any star with ${{\rm{\Sigma }}}_{n}$ ≥ ${{\rm{\Sigma }}}_{\mathrm{cr}}$, its 10 nearest neighbors that also show ${{\rm{\Sigma }}}_{n}$ ≥ ${{\rm{\Sigma }}}_{\mathrm{cr}}$ are its friends; friends of its friends are also friends, and so on. To avoid spurious detections, we require additionally that each group should contain at least ${N}_{\min }^{\star }$ UMS stars. Smaller values of ${N}_{\min }^{\star }$ would lead to more spurious detections, but larger values of ${N}_{\min }^{\star }$ would miss small groups that are not massive enough to contain sufficient number of UMS stars at the high-mass end of their IMF. As we shall see in Section 6, this effect is significant for groups of sizes smaller than 8 pc. In this section, however, this effect does not affect the conclusions reached; thus, we arbitrarily set ${N}_{\min }=10$ in the following analysis. The group identification process is repeated on each level of critical surface density, and we keep track of the connections between "parent" groups on lower levels and "child" groups on higher levels.

The groups found at the first four levels with log(${{\rm{\Sigma }}}_{\mathrm{cr}}\cdot$pc²) = −1.7, −1.5, −1.3, and −1.1 are shown in Figure 4. To help distinguish different groups, each one is enclosed with its convex hull, shown as thick-lined polygons in the figure; 39 groups are found on the lowest level with log(${{\rm{\Sigma }}}_{\mathrm{cr}}\cdot$pc²) = −1.7. Most prominent are the two large groups, labeled "A" and "B" in Figure 4(a), residing in the central part of the tile. Group A consists of the 30 Dor complex, LH 100-LH 99-LH 90, and the diffuse LH 96, LH 111 to its south and east, which have already been mentioned in Section 4. On the other hand, Group B is composed of LH 101–LH 104 and LH 103 in the star-forming regions N158 and N160, respectively. 48 groups are found on the second level with log(${{\rm{\Sigma }}}_{\mathrm{cr}}\cdot$pc²) = −1.5, and the groups are generally smaller in size compared with those on the first level. Group A now splits into two large groups, corresponding to LH 100–-LH 99–LH 90 and LH 96, a smaller LH 111 in the east, and several even smaller ones among them. Group B also splits, with LH 101, LH 104, LH 103, and many other groups gaining their independence. 45 and 40 groups are found on the third and fourth levels, respectively. Results on higher levels are not shown here, but similar processes are continuously repeated, i.e., larger groups vanish and smaller groups emerge.

Zoom In Zoom Out Reset image size

To better display the relations between groups on different levels, we show all the convex hulls of groups found on the first four levels in Figure 5, and the dendrograms of Group A and Group B in Figure 6. One can, of course, construct a single dendrogram containing all groups found on all levels. Here we show only the dendrograms starting from Group A and Group B, which are its "child dendrograms." Both Figures 5 and 6 show that the stellar groups on lower levels very frequently split up into several smaller ones on higher levels. Thus, the stellar structures display a high degree of sub-grouping. This result agrees well with previous works, which also reveal hierarchically subgrouped young stellar structures in individual star-forming regions (e.g., Schmeja et al. 2009; Kirk & Myers 2011; Gouliermis et al. 2014) or over the whole galaxies (e.g., Gouliermis et al. 2010, 2015; Gusev 2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The stellar groups identified in Section 5 have a wide range of sizes, from the largest Groups A and B, which are more than 100 pc in extent, to the smallest ones on parsec scales. Figure 7 shows their cumulative size distribution, which includes stellar groups found on all levels of critical surface densities. The size of a group is estimated using

$$\begin{eqnarray}&&R=\displaystyle \frac{{\displaystyle \sum }_{i}| {x}_{i}-{x}_{c}| {w}_{i}}{{\displaystyle \sum }_{i}{w}_{i}},\end{eqnarray} \tag{ 3 }$$

where x_i is the position vector of the $i\mathrm{th}$ group member, w_i its weight assigned in the same way as described in Section 3 to account for the local photometric completeness, and x_c the group center, which is calculated as

$$\begin{eqnarray}&&{x}_{c}=\displaystyle \frac{{\displaystyle \sum }_{i}{x}_{i}{w}_{i}}{{\displaystyle \sum }_{i}{w}_{i}}.\end{eqnarray} \tag{ 4 }$$

The cumulative size distribution is approximately a power law between ∼10 pc to >100 pc, with significant flattening at several parsecs. As mentioned in Section 5, stellar groups are identified on levels of critical surface densities spaced by equal logarithmic steps of 0.2 dex; with smaller steps and thus more levels, more stellar groups can be included in obtaining the cumulative size distribution. We carried out a test by reducing this step to 0.1 and 0.05 dex and repeating the group identification process; the shapes of the resultant cumulative size distributions, and specifically the slopes of their power-law parts, are almost unchanged (not shown), suggesting that this does not affect the conclusion reached in this section. On the other hand, the flattening is due to incompleteness in our group-identifying algorithm. As mentioned in the previous section, smaller groups may be missed if they contain fewer UMS stars than the required minimum number of stars, ${N}_{\min }^{\star }$. Thus, we have calculated the cumulative size distribution multiple times for stellar groups identified with ${N}_{\min }^{\star }=5$, 10, and 15. As expected, fewer small-sized groups are identified with increasing values of ${N}_{\min }^{\star }$. Beyond R = 8 pc, however, the cumulative size distribution remains almost unaltered. Thus, we suggest that the stellar groups found by our group-identifying algorithm are complete above R = 8 pc.

Zoom In Zoom Out Reset image size

In the range not affected by incompleteness, the cumulative size distribution can be well fitted with a single power law with power-law index α = −1.6 ± 0.3 (fitting error). This indicates a scale-free behavior and agrees with the scenario of hierarchical star formation (Elmegreen et al. 2000). On the other hand, the cumulative size distribution for substructures inside a fractal follows

$$\begin{eqnarray}&&N(\gt R)\propto {R}^{-D},\end{eqnarray} \tag{ 5 }$$

where N is the number of substructures with sizes larger than R, and D is the fractal dimension (Mandelbrot 1983; Elmegreen & Falgarone 1996). Thus, the UMS stars have a fractal dimension of D₂ = 1.6 ± 0.3, by comparing the power-law cumulative group size distribution with Equation (5). Note that D₂ is the fractal dimension of the UMS stars' ${projections}$ on the two-dimensional plane perpendicular to the line of sight. It is not easy to recover the fractal dimension of the UMS stars' ${volume}$ distributions, D₃, from D₂. Beech (1992) suggested a relation D₃ = D₂ + 1; however, this relation holds only if the perimeter-area dimension of a projected three-dimensional structure is the same as that of a slice (Elmegreen & Scalo 2004), and moreover, this relation is often shown to be unrealistic for star-forming regions (e.g., Sánchez et al. 2005). Based on simulations, Gouliermis et al. (2014) proposed a more general conversion between D₃ and D₂, which would suggest D₃ = D₂ = 1.6 for this case; but this may still remain uncertain, since this real star-forming complex can be different in its three-dimensional structures from, and maybe also more complicated than, their simulated fractals.

The ISM also exhibits hierarchical structures of clumps and filaments on all scales (Larson 1981; Elmegreen & Scalo 2004). These structures have power-law cumulative size distributions (e.g., Elmegreen & Falgarone 1996; Roman-Duval et al. 2010), very similar to the UMS stars as found here. Typical fractal dimensions reported are close to D₃ = 2.4; if D₂ = ${D}_{3}-1$  = 1.4 holds for the ISM, this value is consistent with that of the UMS stars, within the errors. This (at least qualitatively) supports a scenario of hierarchical star formation, where the newly born stars follow the the gas distribution (Elmegreen et al. 2000). For the origin of hierarchical structures in the ISM, the role of turbulence as the driving source is highlighted by recent studies (e.g., Federrath et al. 2009); as suggested by Elmegreen et al. (2000), however, this is also possibly accompanied by agglomeration with fragmentation (Carlberg & Pudritz 1990; McLaughlin & Pudritz 1996) as well as self-gravity (de Vega et al. 1996).

Beyond the regime of binary and multiple systems on very small scales, the young stars in a few Galactic star-forming regions are hierarchically clustered with fractal dimensions D₂ = 1.4, as derived by Larson (1995) for Taurus (see also Gomez et al. 1993), or D₂ = 1.5, as estimated by Simon (1997) for Taurus, Orion, and Ophiuchus. The projected fractal dimension we have derived agrees well with these studies. Gouliermis et al. (2014) reported bimodal stellar clustering for the SMC star-forming region, NGC 346; the auto-correlation function (ACF) of young stars is well described by a broken power law, with a shallower slope for scales smaller than ∼5.8 pc and a steeper one for scales above this value. The slope of the latter part corresponds a projected fractal dimension of D₂ = 1.4, also consistent with our result. They have shown that the ACF is best explained by two distinct stellar clustering components, with a centrally concentrated one dominant in the central part of the region, and an extended, hierarchical distribution across the observed field. In the case of 30 Dor-N158–N159–N160, however, while there may be centrally concentrated, non-hierarchical components on small scales (e.g., star clusters), their existence is not revealed by the cumulative size distribution shown here.

Significant for the intense star formation, the 30 Dor-N158–N159–N160 complex has been proposed to arise from galaxy interactions. de Boer et al. (1998) suggested that the star formation in this area may have been induced by the bow shock formed as the LMC moves through the Milky Way's halo. Alternatively, Fujimoto & Noguchi (1990) and Bekki & Chiba (2007) considered the interaction between the Magellanic Clouds. Despite these possible external influences, however, no significant difference in the value of D₂ is found compared with the star-forming regions mentioned above, given the measurement uncertainties.

On larger scales, young stars also exhibit galaxy-wide hierarchical distributions. Specifically, Elmegreen & Elmegreen (2001), Elmegreen et al. (2006, 2014), and Gouliermis et al. (2015) reported power-law cumulative size distributions of the young stellar groups in a sample of galaxies; in contrast, Bastian et al. (2007, 2009) suggested that the young stellar groups in M33 and the LMC have log-normal size distributions. However, if we deem the flattening on small scales in their cumulative size distributions (their Figures 5 and 1, respectively) as caused by incompleteness, the large-scale parts can also be well described by single power laws. Note that Bastian et al. (2007, 2009) required each of their identified stellar groups should contain at least five sources, which may miss small-sized groups; but they did not discuss the issue of incompleteness as has been done earlier in this section.

On the other hand, the reported projected fractal dimensions of the galaxy-wide hierarchies span a range from close to 1.0 to as large as 1.8 or more (Elmegreen & Elmegreen 2001; Elmegreen et al. 2006, 2014; Gouliermis et al. 2015). D₂ derived in this work is close to those for NGC 628 (D₂ = 1.5; Elmegreen et al. 2006), NGC 6503 (D₂ = 1.7; Gouliermis et al. 2015), and part of the galaxy sample of Elmegreen et al. (2014), but it can deviate more than the measurement uncertainties compared with other ones. Elmegreen et al. (2014) found that, compared with large spiral galaxies or low surface brightness dwarfs, the starburst dwarfs or HII galaxies have larger projected fractal dimensions if they contain one or two dominant young stellar complexes; thus, the difference in D₂ may reflect different clustering under different conditions. Second, older stellar populations have a smaller degree of substructures due to dynamical evolutions (Gieles et al. 2008; Bastian et al. 2009; Parker et al. 2014; Gouliermis et al. 2015); this effect may also contribute to the difference in D₂, since these studies rely on different stellar samples that may be complicated mixtures of populations of various ages. Moreover, it is worth noting that the derived fractal dimensions may also depend on the adopted method (Federrath et al. 2009), but in this work we do not attempt to explore this effect.

Using the same VMC tile as in this work, Romita et al. (2016) reported the discovery of abundant embedded star clusters in molecular clouds. The size distribution of the embedded clusters spans the range of 0.25–2.25 pc with a peak at ∼1.1 pc, very different from ours. Through their visual inspections, they preferentially identified small and compact cluster candidates, which have high contrast of stellar surface density with respect to their surrounding areas. The group identification algorithm used in this work, however, does not introduce a preferential scale and reveals star groups over a continuous range of sizes. On the other hand, as mentioned above, the star groups identified with the relatively rare UMS stars are incomplete below R = 8 pc. Thus, Romita et al. (2016) and this work focus on different subsections of the full hierarchy of young stellar structures. Bonatto & Bica (2010) investigated the size distributions of star clusters and "non-clusters," which are essentially nebular complexes and stellar associations. The former class has a steeper distribution, and the latter has a shallower distribution than ours. Their results are different from ours owing to two reasons. First, they have treated the two classes separately, while we do not try to distinguish star clusters, associations, or larger stellar structures. Second, they have used the catalog of extended objects of Bica et al. (2008), which, in turn, was compiled from many previous works. Thus the effect of selection bias is hard to assess.

In this work we investigate the hierarchical stellar structures of the 30 Dor-N158–N159–N160 star-forming complex traced by young UMS stars observed with the VMC survey. We calculate the surface densities of the UMS stars, whose values cover two orders of magnitudes. We identify groups of UMS stars on different levels of critical surface densities. The larger-sized groups on lower levels often contain several smaller-sized ones on higher levels, and a high degree of sub-grouping is revealed by the dendrograms constructed to show the "parent–child" relations between the groups on different levels. The stellar groups have sizes ranging continuously from several parsecs to more than 100 pc. Without characteristic sizes, the stellar groups have a power-law cumulative size distribution, with the power-law index indicating a projected fractal dimension D₂ = 1.6 ± 0.3. We suggest that the cumulative size distribution of the UMS stars is related to their parental ISM in a scenario of hierarchical star formation, in which newly born stars follow the gas distribution. The projected fractal dimension derived here is close to those reported for other star-forming regions, as well as for several galaxies; the difference with respect to some other galaxies is also discussed.

We are grateful to the anonymous referee for the detailed and constructive comments on this paper. The analysis in this article is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere (ESO) under ESO program 179.B-2003. We thank the Cambridge Astronomical Survey Unit (CASU) and the Wide Field Astronomy Unit (WFAU) for providing calibrated data products under the support of the Science and Technology Facility Council (STFC) in the UK. We thank Jim P. Emerson for helpful suggestions. N.-C.S. and R.d.G. acknowledge funding support from the National Natural Science Foundation of China through grants 11373010, 11633005, and U1631102.