The VMC Survey. XXIX. Turbulence-Controlled Hierarchical Star Formation in the Small Magellanic Cloud

In this paper we report a clustering analysis of upper main-sequence stars in the Small Magellanic Cloud, using data from the VMC survey (the VISTA near-infrared YJKs survey of the Magellanic system). Young stellar structures are identified as surface overdensities on a range of significance levels. They are found to be organized in a hierarchical pattern, such that larger structures at lower significance levels contain smaller ones at higher significance levels. They have very irregular morphologies, with a perimeter-area dimension of 1.44 +/- 0.02 for their projected boundaries. They have a power-law mass-size relation, power-law size/mass distributions, and a lognormal surface density distribution. We derive a projected fractal dimension of 1.48 +/- 0.03 from the mass-size relation, or of 1.4 +/- 0.1 from the size distribution, reflecting significant lumpiness of the young stellar structures. These properties are remarkably similar to those of a turbulent interstellar medium (ISM), supporting a scenario of hierarchical star formation regulated by supersonic turbulence.


INTRODUCTION
Star formation leads to hierarchical young stellar structures ( 10 8 yr), including star clusters, associations, and complexes of increasing size and decreasing density (Ivanov et al. 1992;Efremov 1995;Bastian et al. 2006. A complex may be subclustered into star clusters and associations (Eigenson & Yatsyk 1988;Piskunov et al. 2006;de la Fuente Marcos & de la Fuente Marcos 2008;Elmegreen 2011), and the associations may in turn contain a number of subgroups and star clusters (e.g. the Orion OB1 association; Genzel & Stutzki 1989;Kuhn et al. 2014;Da Rio et al. 2014).
The young stellar structures are a natural outcome of star formation, which itself is a hierarchical process and appears correlated in both space and time. This is supported by the age difference-separation relation of star clusters Grasha et al. 2017a), young star clusters forming in pairs or groups (Bhatia & Hatzidimitriou 1988;Hatzidimitriou & Bhatia 1990), as well as the positional correlations of newly-formed stars (Gomez et al. 1993;Larson 1995;Gieles et al. 2008;Gouliermis et al. 2015) and star clusters (Zhang et al. 2001;Scheepmaker et al. 2009;Grasha et al. 2017b). In-terstellar turbulence may play a key role in regulating hierarchical star formation across a wide scale range , possibly aided by fragmentation with agglomeration and self gravity (Carlberg & Pudritz 1990;de Vega et al. 1996). The young stellar structures undergo rapid dynamical evolution after their formation. Most of the newly formed star clusters are short-lived, and only a small fraction are able to survive the early disruptive processes (Lada & Lada 2003;Portegies Zwart et al. 2010). On larger scales, stellar associations and complexes are not gravitationally bound, and they dissolve into the general galactic field within no more than several hundred million years (Efremov 1995;Gieles et al. 2008;Gouliermis et al. 2015). It is therefore important to study their properties, formation and evolution, which will provide insights into a galaxy's star formation and dynamical processes on multiple scales.
The star cluster population of the SMC has been extensively investigated (e.g. Pietrzyński et al. 1998;de Grijs & Goodwin 2008;Gieles et al. 2008;Piatti et al. 2016). In contrast, its stellar associations have received far less attention [see Table 4 of Bica et al. (2008) for a catalog of SMC associations]. Note that the SMC is too small to contain many stellar complexes. The young stellar structures in the SMC are thus poorly understood especially in terms of their hierarchical nature.
The early work of Hodge (1985) reported 70 OB associations in the SMC. With a mean diameter of 77 pc, their SMC associations have a similar mean size as those in the LMC, but are smaller than those in the Milky Way, M31, and M33. As pointed out by the author, however, this may not be a true size difference and might result from selection effects -the different abilities to recognize associations in galaxies at different distances, and the different criteria, often subtle and subjective, to separate the stellar distribution into associations. Battinelli (1991) proposed an objective method, which they used to identify associations in the SMC. They reported a mean diameter of 90 pc for the associations, which is comparable to that of Hodge (1985). However, as Bastian et al. (2007) pointed out, this characteristic size may also be a selection bias: the mean size would change when adopting different values for the breaking scale and minimum number of stars in the identification algorithm. Actually, the stellar associations may not have any characteristic sizes. Gieles et al. (2008) investigated the spatial distributions of stars of various ages in the SMC. They found that stars are born with significant substructures in their spatial distribution with a fractal dimension of 2.4. Older stars become successively less subclustered, and the stellar distribution becomes smooth on a timescale of ∼75 Myr. Unfortunately, they did not investigate the properties of the young stellar structures in greater detail. Bonatto & Bica (2010) studied the size distributions of star clusters and "non-clusters" (nebular complexes and stellar associations) in the Magellanic Clouds. They found that the size distributions follow power laws, which is in agreement with a scenario of hierarchical star formation. However, they relied on the catalog of extended sources of Bica et al. (2008), which is in turn based on the extensive previous catalogs from different authors. Thus, it is not easy to assess the effects of selection biases and incompleteness.
It is therefore not trivial to identify young stellar structures 1 while avoiding subjectivity and selection bias. One important reason comes from their hierarchical nature, in the sense that they have irregular morphologies and abundant substructures -they themselves may be the substructures of larger ones as well (see also the discussion by ). This calls for 1 In the following we shall not distinguish star clusters, associations, and complexes, and we choose to refer to all of them as young stellar structures. systematic, objective and unbiased identification algorithms. A contour-based map analysis technique, proposed by Gouliermis et al. (2015Gouliermis et al. ( , 2017, has proved to be effective to achieve this end. The idea is to identify young stellar structures as overdensities on a series of significance levels from the surface density map of young stars. This method provides an intuitive view of the young stellar structures and, in particular, their hierarchical nature. This technique, or similar ones based on the same principles, have been applied to a number of nearby galaxies (e.g. Gouliermis et al. 2010Gouliermis et al. , 2015Gouliermis et al. , 2017 and star-forming complexes (e.g. Sun et al. 2017a,b). Still, it has not yet been applied to the SMC as a whole.
Considering that the young stellar structures span multiple scales, the observations require both high spatial resolution and wide-field coverage. The VMC survey (Cioni et al. 2011) provides a precious opportunity to meet this challenge. With a spatial resolution of 1 ′′ (0.3 pc at the SMC's distance), it covers both Clouds' full classical extents at B = 25 mag arcsec −2 (Bothun & Thompson 1988), the Bridge between the two galaxies, with two additional fields centered on their associated Stream. Working at near-infrared wavelengths, the VMC survey is less susceptible to the effects of reddening, which is usually important for the young stellar populations. In this paper, we apply the contour-based map analysis to the SMC, using upper main-sequence (UMS) stars from the VMC survey. Our goal is to understand the hierarchical properties of its young stellar structures, which will also demonstrate their formation mechanism(s). In Section 2 we outline the VMC data used in this work, and Section 3 describes how we build the stellar density map. Young stellar structures are identified in Section 4. We analyze their properties in Section 5, followed by a discussion in Section 6. We finally close this paper with a summary and conclusion, leaving the issue of their evolution for a future paper.

SELECTED VMC DATA
Data used in this work come from the VMC survey, which is carried out with the VISTA telescope (Sutherland et al. 2015). We refer to Cioni et al. (2011) for a comprehensive description of the VMC survey. The original pawprints 2 are processed by the VISTA Data Flow System (VDFS; Irwin et al. 2004), and PSF photometry is performed based on stacked tiles combining all good-quality epochs (Rubele et al. 2012(Rubele et al. , 2015. We use PSF photometry in this work, since it suffers less from source crowding than aperture photometry (Tatton et al. 2013). The data are retrieved from the VISTA Science Archive 3 (VSA; Cross et al. 2012).
There are 27 SMC tiles in the VMC survey. In this paper, we only use the central 4 × 4 tiles covering the main body of the SMC and its close vicinity. The other tiles are in the more outer regions and have very low stellar surface densities; thus they are not analyzed in this work. The names and extents of the 16 tiles for analysis are shown in Fig. 1 overlaid on a stellar density map.
2 Because there are gaps between the detectors, VISTA observes a contiguous area of sky using a sequence of six offset observations. Each of the six observations is called a pawprint, while the combined image, covering a contiguous area of ∼1.5 deg 2 , is referred to as a tile. Tiles SMC 3 3, 3 5, 4 3, 5 2, and 5 4 are part of Data Release 4 while the other tiles are currently still proprietary to the VMC team (they will be included in Data Release 5 expected in 2018). The stellar density map of Fig. 1 is obtained from simple star counts of VMC sources with K s < 20 mag in 0.01 • × 0.01 • cells. This magnitude cut includes the UMS, red giant branch (RGB), and the red clump (RC) of the SMC. The SMC's main body is clearly revealed in the central region of the map. Tiles SMC 3 5, 4 5, and 5 5 cover the inner part of the SMC's Wing, which is located to the east of the main body. It is obvious that the 16 tiles contain the bulk of the SMC's stars.
Adjacent tiles partly overlap with each other; as a result, there are duplicate sources in the combined source catalog. If a source is located in the overlap region and has detections in two or more tiles, we keep only the detection from the tile whose center is closest to the source's position. In this way, the duplicate sources are removed and all sources are unique in the combined catalog. According to the tiling algorithm of VISTA, strips along the upper and lower outer edges of each tile are not covered to the same depth as the bulk of the tile. This is because only one pawprint contributes there, while the bulk of the tile is covered by 2-6 pawprints. Ideally, one needs to coadd the objects in the overlap regions between two tiles to obtain uniform depth. However, this effect does not make any significant difference in the context of this work, since the sources under investigation are very bright compared with the survey's magnitude limits (see Section 3.1).
Two overdensities are located in tiles SMC 5 2 and 6 4. They correspond to the Galactic globular clusters, 47 Tucanae (47 Tuc) and NGC 362. In the following analysis, we exclude all sources within 0.5 • and 0.2 • from their centers 4 , respectively, which can effectively eliminate the influence of the two globular clusters. The stellar density map also shows a vertical strip devoid of stars at ∆R.A. ∼ 0.3 • . This arises from a gap between tiles SMC 5 3 and 5 4. Since SMC 5 3 is slighted offset to the west, this gap is not covered by the observations. A complementary VISTA observation is ongoing to fill this gap. In this work, we will use the Magellanic Cloud Photometric Survey (MCPS; Zaritsky et al. 2000Zaritsky et al. , 2002 to fill this gap, as detailed in Section 3.3. Figure 2 shows the (J − K s , K s ) color-magnitude diagram (CMD) of all sources in the 16 selected SMC tiles. The main sequence (MS), RGB, and RC are all clearly visible. Brighter than the tip of the RGB is the asymptotic giant branch (AGB), and the red supergiant branch (RSG) is located slightly bluer than the RGB and AGB. Redward of J − K s = 1.0 mag are background galaxies (BGGs), and the vertical spurs at J − K s = 0.75 mag and 0.35 mag are foreground stars (FGSs) in the Milky Way.

Selection of UMS Stars
In the right-hand panel of Fig. 2, we also show the PARSEC (version 1.2S; Bressan et al. 2012) isochrones of ages log(τ /yr) = 7.0, 8.0, and 9.0. The solid and dashed isochrones correspond to 20% and 10% of solar metallicity, respectively 5 . The former metallicity is derived from H II regions and young stars (Russell & Dopita 1992), while the latter is typical of the RGB population in the SMC (Dobbie et al. 2014). Offsets of 0.026 mag in J and 0.003 mag in K s have been subtracted from the isochrones to convert the model magnitudes from the Vega system to the VISTA system (Rubele et al. 2012(Rubele et al. , 2015. The isochrones are shifted by an SMC distance modulus of (m − M ) = 18.96 mag (de Grijs & Bono 2015) and corrected for a foreground Galactic extinction of A V = 0.12 mag (Schlegel et al. 1998). We adopt the extinction coefficients A J /A V = 0.283 and A Ks /A V = 0.114, which are computed from the Cardelli et al. (1989) extinction curve with R V = 3.1 (Girardi et al. 2008). The extinction coefficients in J and K s depend only very weakly on the stellar spectral types or the adopted extinction laws (see e.g. van Loon et al. 2003).
UMS stars are selected based on 13.0 < K s < 18.0 mag and −0.3 < J − K s < 0.2 mag, which is indicated by the black solid box in Fig. 2. The isochrones suggest that the selected stars are younger than 1 Gyr. Since the VMC survey has typical saturation limits of J = 12.7 mag and K s = 11.4 mag (Cioni et al. 2011), the upper limit of K s = 13.0 mag ensures that the selected stars are not saturated in either band. This limit corresponds to a stellar mass of ∼19 M ⊙ for an age of 10 Myr and a 20% solar metallicity. The color range is chosen with consideration of interstellar extinction. In addition to the foreground Galactic extinction, the stars also suffer from extinction inside the SMC. This component is, however, very difficult to model theoretically or derive observationally. First, extinction in star-forming regions often exhibits 5 The PARSEC isochrones have adopted a solar metallicity of Z ⊙ = 0.0152, and are computed following a Y = 0.2485 + 1.78Z relation for the helium abundance. Here Y and Z are the mass fractions of helium and metal elements. . In the right-hand panel, PARSEC (version 1.2S) isochrones are overplotted, shifted by a distance modulus of (m − M ) 0 = 18.96 mag and a foreground Galactic extinction of A V = 0.12 mag. They correspond to ages of log(τ /yr) = 7.0 (black), 8.0 (red), and 9.0 dex (blue) with 20% or 10% solar metallicity (solid or dashed lines, respectively). The black box shows the criteria adopted for selecting UMS stars. The color scales in both panels are identical; the labels, isochrones, and selection box are shown separately simply for clarity. significant spatial variations (e.g. Lombardi et al. 2010). Second, dust may lie in the foreground or background of a star, and only the foreground dust contributes to its extinction. Thus, an extinction map derived from a group of stars is not always directly applicable to another group of stars, if they have different distributions along the line of sight with respect to the dust. Moreover, extinction in the Magellanic Clouds exhibits a population dependence (Zaritsky 1999). Extinction has a very small effect on the infrared magnitudes, but may shift some UMS stars redward out of the selection box. Thus, we have adopted a wide color range in the selection criteria to include the reddened stars as completely as possible. For instance, an extinction of A V = 1.8 mag is needed to shift a star redward from the log(τ /yr) = 7.0 dex isochrone out of the selection box; such a high extinction is very rare in the SMC (see e.g. Haschke et al. 2011;Rubele et al. 2015). The UMS sample contains 46,148 stars.
Because we have adopted a wide color range in the selection, the UMS sample consists of both UMS stars and slightly evolved stars that are moving redward off the MS or on the blue loops. However, their number is likely small, since a star does not spend much time on these rapid evolutionary phases. On the other hand, many of these stars may still have very young ages; if so, they can also help trace the young stellar structures. Rubele et al. (2018) analyzed the star formation history (SFH) of the SMC based on a stellar population synthesis technique. Their mapped region is the same as that adopted in this work. They found that the SMC had a very low star-forming rate at log(τ /yr) = 8.5-9.0 (lookback time), while more recently it has experienced increased star-forming activity. Based on their SFH, it is possible to assess the age distribution of stellar populations in the SMC (see their Section 5.4). We have done such an analysis for the UMS sample based on their results. We found that the sample has a median age close to log(τ /yr) = 7.7 (50 Myr); approximately 72%, 86%, and 96% of the selected stars are younger than log(τ /yr) = 8.0, 8.2, and 8.4 (100, 160, and 250 Myr), respectively. Thus, while the comparison with isochrones ( Fig. 2) gives an upper age limit of 1 Gyr, a detailed SFH analysis suggests a much younger age distribution for the UMS sample.
Here we also make an important distinction between the ages of the selected UMS stars and the ages of the young stellar structures. In this work, the young stellar structures will be identified as surface overdensities (Section 4). On the other hand, recent papers have shown that young stars are the most subclustered while old ones have very smooth spatial distributions (Gieles et al. 2008;Gouliermis et al. 2010;Sun et al. 2017b). This is because the young stellar structures undergo rapid, disruptive dynamical evolution and/or because of the spatial overlap of different generations of star formation (Elmegreen 2018). As a result, the surface overdensities are dominated by young rather than old stars. Gieles et al. (2008) analyzed the age-dependent twopoint correlation functions of the SMC's overall stellar distributions. They found that the stellar distribution becomes indistinguishable from a uniform distribution at an age of 75 Myr. Although this value cannot be taken as an upper age limit for a specific individual stellar structure, statistically speaking most stellar structures should be dispersed roughly on this timescale (oth-erwise, the stellar distribution would still be highly nonuniform).   On the other hand, most young stellar structures are unbound, except for some compact star clusters on small scales. They could disperse due to the stellar random motions, which are related to the turbulent motions in their natal star-forming clouds. The turbulent crossing time is typically several tens of millions of years for a cloud of 100 pc (see e.g. . The structure dispersion timescale should be comparable to or a few times the turbulent crossing time, if stellar random motions dominate their dispersion process (this is a valid assumption especially for large-scale, loosely bound or unbound structures). In addition, the recent LMC-SMC encounter at 200-400 Myr ago must have had a strong impact on the SMC (Maragoudaki et al. 2001;Bekki & Chiba 2007); some of the stellar structures older than this age could have been more easily destroyed under the LMC's tidal force. Considering all this, it is reasonable to assume that most young stellar structures have an upper age limit of the order of 100 Myr.
One may still question whether the results can be affected by contamination of stars in the UMS sample which are older than the typical ages of young stellar structures. We have done a test with two additional UMS samples, each selected with a different fainter magnitude cut at K s = 17.5 or 18.5 mag. They have different stellar age distributions compared with the original sample. 77% or 58% of their stars are younger than 100 Myr, respectively, so they have a reduced or increased fraction of older stars. We then repeat the entire analysis as described below; the results do not change significantly and the derived quantities (perimeter-area dimensions, slopes of size/mass distributions, fractal dimension, etc.; see Section 5) are consistent within the uncertainties. This test again confirms that the young stellar structures are dominantly revealed by the young rather than the old stars. Thus, the conclusions reached below are robust and not affected by the contamination by old stars.
In the sample selection we have adopted a single distance modulus for the SMC. However, the SMC is known to be significantly elongated by more than 25-30 kpc along the line of sight (e.g. Subramanian & Subramaniam 2009;Ripepi et al. 2017). As will be shown later, our identified young stellar structures lie primarily in the central main body of the SMC. The SMC's line-of-sight depth is much smaller in this limited area. On the other hand, Ripepi et al. (2017) showed that the young classical Cepheids (< 140 Myr) are distributed in a geometry that is less elongated than the older ones. Inspection of their Fig. 15 suggests that most of the young Cepheids are located at 60-67.5 kpc in the SMC's central area. This corresponds to a −0.07/+0.19 mag difference from the adopted distance modulus. The potential zero-point difference is small with respect to the size of the UMS box; thus, we do not consider the effect of SMC's line-of-sight depth in our analysis.

The KDE Map
We construct a surface density map of the selected UMS stars using kernel density estimation (KDE). This is done by convolving the stellar distribution with a Gaussian kernel. In principle, the choice of the kernel width is arbitrary and depends on the scale of the object we are interested in. A narrow kernel preserves the small-scale structures in the map but may suffer from large statistical fluctuations. On the other hand, a wide kernel reduces statistical noise, but the resultant KDE map would have poor spatial resolution. In this section, we try to build a KDE map where the small-scale structures are still preserved. As a result, we look for a kernel that is as small as possible. Based on tests with different widths, we found that a Gaussian kernel with a width of 10 pc (standard deviation) is appropriate to balance spatial resolution and noise. With this kernel, it is possible to resolve structures down to scales of ∼10 pc. The corresponding KDE map is shown in Fig. 3.
A number of well-known features are evident in the KDE map. The majority of UMS stars are distributed in the main body of the SMC, which is northeast-southwest elongated and also known as the Bar. A diffuse Wing of stars resides to the east of the Bar, and another stellar extension is seen in the extreme SW (SW Extension). To the northeast there is a diffuse Shell of stars perpendicular to the Bar. All of these features are labeled in the figure, and they agree well with previous studies by, e.g., Cioni et al. (2000b) and Zaritsky et al. (2000).
There is a vertical strip devoid of stars at ∆R.A. ∼ 0.3 • , which arises from the gap between VMC tiles SMC 5 3 and 5 4 (see Fig. 1). However, our map analysis requires a contiguous spatial coverage. Thus, we use data from the MCPS survey (Zaritsky et al. 2002(Zaritsky et al. , 2004 to fill this gap, which is detailed in the following subsection.

Filling the Gap with MCPS
The MCPS survey is a photometric survey of the Magellanic Clouds in the U , B, V , and I bands. The limiting magnitude is approximately V = 20 and I = 21 mag, depending on the local crowding of the images. To fill the gap, we first define a rectangular region of 0.45 • × 1.40 • centered on the gap. The region is shown as the solid rectangle in Fig. 3 and is referred to as Region A hereafter.
Next, we select a sample of UMS stars in Region A from the MCPS data. Figure 4 shows the (B − V , V ) CMD of the MCPS sources in this region. The MS, RGB, and RC are all clearly seen. We cross-matched the MCPS sources with the UMS stars selected in Section 3.1, using a matching radius of 1 ′′ . The matched sources are shown as the contours in Fig. 4. Thus, we use the red polygon, which approximately encloses the contours, to select a sample of UMS stars based on MCPS photometry. This "MCPS UMS sample" contains 13,991 stars in total. Note that extinction in this region is relatively low and is not a matter of concern here (see e.g. Fig. 4 of Haschke et al. 2011).
Third, a KDE map of Region A is constructed for the MCPS UMS sample in the same way as in the previous subsection. Note that the pixel values of the VMC and MCPS maps are not necessarily the same since they are based on different samples. Figure 5 shows the pixel-topixel comparison of the VMC and MCPS maps of Region A. The comparison is done only with the "good" pixels -pixels in the gap, within 10 pc (or 34 ′′ .4, the kernel width of the KDE maps) from the gap edge, or within 10 pc from the region edge are excluded from the comparison. Figure 5 shows that there is a tight correlation between the pixel values of the VMC and MCPS maps. The correlation can be fitted with a linear function (solid line in Fig. 5),  Fig. 6. In general, they agree with each other very well, except for areas affected by the gap or by the region edge. This ensures that we can fill the gap with MCPS data. Finally, we try to combine the two maps. A weight map is constructed, where the pixels have a weight of w = 1.0 if they are more than 30 pc (three kernel widths) from the gap edge; a weight of w = 0.0 is assigned to pixels inside the gap or within 10 pc from the gap edge; linear interpolation is used to assign weights to all other pixels. The VMC and MCPS maps are then combined as where µ comb is the pixel value of the combined map. Adopting this method, the combined map (Fig. 6, righthand panel), from the outer areas to the gap, changes smoothly from the VMC values to the (calibrated) MCPS values, effectively filling the gap. The complete gap-filled map is shown in Fig. 7. The map has a median value of 2 × 10 −5 pc −2 , a mean value of 0.002 pc −2 , and a standard deviation of 0.004 pc −2 . Our identification and analysis of the young stellar structures is based on this map.

IDENTIFICATION OF YOUNG STELLAR STRUCTURES
We use the contour-based map analysis technique proposed by Gouliermis et al. (2015Gouliermis et al. ( , 2017 to identify young stellar structures in the SMC. The idea is to identify them as surface overdensities over a range of significance levels. In the KDE map (Fig. 7), the iso-density contours are obtained from 1σ to 15σ in equal steps of 1σ (σ = 0.004 pc −2 ; note that the median value of the KDE map is very close to zero compared with the standard deviation). On each significance level, any iso-density contour enclosing an overdensity is regarded as a can-  didate young stellar structure. The iso-density contours are thus regarded as the (projected) "boundaries" of the candidate young stellar structures. We next determine physical parameters of the candidates. The size of a structure (R) is estimated with the radius based on a circle of the same area as that of its boundary. The structure's mass is characterized by the number of UMS stars inside its boundary (N * ). To first-order approximation, N * is proportional to the mass of the structure, assuming the structures have similar ages and that the stellar initial mass function (IMF) is fully sampled (Gouliermis et al. 2015;Sun et al. 2017b; see also Section 5.3 for a more detailed discussion). N * is obtained from simple star counts of the VMC UMS sample. For the gap of VMC observations, we use instead the MCPS UMS sample (according to Eq. 1, each MCPS-selected UMS star accounts for 0.36 of a VMCselected one). We will hereafter refer to N * as the "stellar number" for simplicity. Finally, the structure's surface density is calculated as Σ = N * /(πR 2 ).
However, the candidates may include spurious detections as well. A commonly adopted method to reject unreliable detections is to require a minimum stellar number (N min ) to define a young stellar structure (e.g. Bastian et al. 2007Gouliermis et al. 2015Gouliermis et al. , 2017Sun et al. 2017a,b). In other words, a candidate is regarded a reliable young stellar structure only if N * ≥ N min . The choice of N min is arbitrary -a large N min can reject most unreliable detections but may remove genuine structures as well; a small N min , on the other hand, leads to more unreliable structures. In the following analysis, we adopt N min = 5 stars, which has been used extensively in previous works. Furthermore, candidates on the lowest significance levels may arise from random density fluctuations. However, we cannot rule out the possibility that they may be genuine structures with low surface densities. In order to reject the less probable candidates of this kind, we regard candidates on the 1σ and 2σ significance levels as genuine structures if their boundaries enclose one or more structures on the 3σ-15σ significance levels.
After rejecting the less reliable candidates, we finally identify 556 young stellar structures. They include 14,30,119,90,72,59,47,45,36,19,14,10,5,3, and 3 structures on the 1σ-15σ significance levels, respectively. The structures are listed in Table 1, along with their IDs, significance levels, coordinates, and physical parameters. In the following, we will refer to a specific structure using its ID number with a prefix "S" (short for "structure"; for instance, "S1" refers to the first structure in Table 1). The boundaries of all identified structures are shown in Fig. 8. Most of the identified structures belong to the Bar region, while another group of structures are from the Wing. Very few structures are identified in the Shell or in the SW extension, since the stellar surface density is very low there. Also note that some (but not all) structures on low significance levels correspond to significant portions of the whole SMC. In fact, the largest structure S1 covers essentially the whole main body of the galaxy. Only on higher significance levels do they break down to independent, internal structures on sub-galactic scales. This is very different from disk galaxies, where the young stellar structures are clearly distinct even on low significance levels (e.g. Gouliermis et al. 2015Gouliermis et al. , 2017. This is possibly because disk galaxies have a more ordered structure, while the SMC is more irregular and its stellar populations are more mixed. Note.

The Hierarchical Pattern
-Column 1: ID number for each young stellar structure; Column 2: significance level; Column 3: right ascension of the center of each young stellar structure, defined as α = (α min + αmax)/2, where α min and αmax are the minimum and maximum right ascension of the iso-density contour of each structure, respectively; Column 4: same as Column 3 but for the declination; Columns 5-7: size, stellar number, and surface density. Only the first 15 records are shown as an example. The complete catalog is available online. Figure 8 shows that the young stellar structures are organized in a hierarchical manner. Generally, a structure on a low significance level may contain one or more Fig. 8.-Boundaries of all identified young stellar structures, with the colors encoded according to their significance levels. S1, whose dendrogram is shown in Fig. 9, is also labeled. Fig. 9.-Dendrogram of the young stellar structures. The red points indicate the structures on different significance levels, and the black links show their "parent-child" relations. The root corresponds to S1 on the 1σ significance level, and the branches and leaves represent its substructures on higher significance levels.
substructures on the higher significance levels within its boundary. To illustrate this behavior intuitively, we show in Fig. 9 the dendrogram of structure S1 on the 1σ significance level. Dendrograms are structure trees showing the "parent-child" relations of young stellar structures on various significance levels. In S1's dendrogram, the "root" corresponds to S1 itself, while the "branches" and "leaves" represent its substructures on the 2σ-15σ significance levels. It is obvious that S1 branches into a lot of substructures on the 2σ significance level, many of which in turn branches into substructures on the 3σ significance level. This fragmentation-like behavior con-tinues to higher significance levels until the top (15σ) significance level. Thus, the young stars in the SMC exhibit a high degree of hierarchical subclustering, which is typical in many star-forming regions and galaxies (e.g. Gouliermis et al. 2010Gouliermis et al. , 2015Gusev 2014;Sun et al. 2017a,b). Figure 8 also shows that many of the young stellar structures have very irregular boundaries, which are significantly different from smooth curves. This morphological irregularity can be characterized by the perimeter- area relation. In Fig. 10 we show the perimeters (P ) and areas (A) of the boundaries of all identified young stellar structures. In the plot, the blue dashed line shows the perimeter-area relation of geometric circles, which is characterized by a power-law slope of α p = 0.5. This relation "clips" the data points at A < 3 × 10 3 pc 2 by coinciding with their lower limit. Data points at A < 3 × 10 2 pc 2 are highly consistent with this relation. This reflects the resolution effect for the small young stellar structures. Since the KDE map (Fig. 7) is based on a kernel width of 10 pc, structures smaller than or comparable to this size will appear nearly roundish because of the kernel smoothing.

Morphological Irregularity
Beyond A = 3 × 10 3 pc 2 , the data points seem not to be strongly "clipped" by the perimeter-area relation of geometric circles, suggesting that resolution effects are no longer important. Note, according to the definition of structure size (Section 4), this limit corresponds to R = 30 pc, or three times the map resolution 7 . The young stellar structures have a power-law perimeterarea relation beyond this limit. The power-law slope is α p = 0.72 ± 0.01 based on a least-squares fit. The perimeter-area dimension, D p , is defined such that P ∝ A Dp/2 (Falgarone et al. 1991). Thus, the boundaries have a perimeter-area dimension of D p = 1.44 ± 0.02. This is significantly larger than their topological dimension of unity, which quantitively reflects the morphological irregularity of the young stellar structures.

Mass-Size and Surface Density-Size Relations
7 As will be shown in Section 5.4, the identification of young stellar structures is complete beyond 10 pc (the KDE map resolution). However, we caution that the perimeter-area relation between 10-30 pc may still be affected by the resolution effect. This is because a structure can be detected as long as its size is beyond the map resolution (and N * ≥ N min in order not to be rejected), but the shape of its boundary can only be revealed when its size is much larger than the resolution. Otherwise, its shape will still appear roundish due to kernel smoothing. In other words, for a given structure, one needs a lower resolution to detect it, but a higher resolution to resolve its shape.  Figure 11 shows the stellar number-size and surface density-size diagrams of all identified young stellar structures. With a Pearson correlation coefficient of −0.59, the surface densities and sizes of the structures do not show any significant correlations. In contrast, the young stellar structures exhibit a tight power-law correlation between their stellar numbers and sizes, with a Pearson correlation coefficient of 0.90. Using a least-squares fit (solid line) we find that the best-fitting slope is κ = 1.48 ± 0.03.
As mentioned in Section 4, the stellar number of a structure is proportional to its mass, if the stellar IMF is fully sampled and if the structures have similar ages. This approximation will be acceptable if the stellar number-to-mass ratio (ξ) does not change significantly among the young stellar structures. To quantitively assess whether this is valid, we carry out numerical simulations of stellar populations with 20% solar metallicity (typical of young stars in the SMC; Russell & Dopita 1992) and with ages from log(τ /yr) = 6.6 to 8.8 in steps of 0.1. For each population, a total of 20,000 stars are sampled randomly from the Chabrier (2001) lognormal IMF. We then assume that 30% of them are unresolved binary systems (Elson et al. 1998;Li et al. 2013), and assign masses to their secondary stars based on a flat distribution of primary/secondary mass ratios from 0.7 to 1.0 (below this limit the secondary star has little influence on the system's overall brightness). This corresponds to an initial mass of ∼1.55 × 10 4 M ⊙ in total for each stellar population. The J and K s magnitudes for the simulated stars are derived using PARSEC isochrones and shifted to SMC's distance modulus. The UMS stars of each stellar population are selected based on the same criterion as used for the real data (Section 3.1), and the stellar number-to-mass ratios are then calculated as ξ = 1000 × N * /(M/M ⊙ ).
We carry out 100 realizations for each age. The mean values of their stellar number-to-mass ratios are shown as the data points in Fig. 12. Starting from ξ ∼ 3 at log(τ /yr) = 6.6, ξ slightly increases to a maximum of 3.3 at log(τ /yr) = 7.1-7.2, and then declines gradually to 1.7 at log(τ /yr) = 8.0, until it reaches 0.8 at log(τ /yr) = 8.4. Within this age limit, ξ differs by at most a factor of ∼4. ξ drops sharply at log(τ /yr) = 8.5 by an order of magnitude and continues to decrease toward older ages. Stellar populations older than log(τ /yr) = 8.8 have no or negligible numbers of stars that meet the UMS selection criteria adopted in this work.
As discussed in Section 3.1, statistically speaking most SMC young stellar structures should have ages younger than or comparable to ∼100 Myr. For these structures, the stellar number-to-mass ratio is almost constant around 2, with minimum and maximum values of 1.7 and 3.3, respectively. Structures with ages log(τ /yr) = 8.0-8.4 (∼100-250 Myr) have lower stellar number-to-mass ratios, but they are still not very different. In Section 3.1 we have mentioned that 96% of stars in the UMS sample are younger than 250 Myr. Thus, we expect very few structures, if any, older than this age. Based on this analysis, we suggest that the variation of ξ caused by evolutionary effects should be small for most SMC young stellar structures.
On the other hand, the error bars in Fig. 12 show the standard deviations of ξ of the 100 realizations for each age. They reflect the variation of ξ caused by the stochas- 1.60 ± 0.08 10 1.7 ± 0.1 11 1.7 ± 0.2 12 1.6 ± 0.2 ≥13 1.6 ± 0.2 total 1.48 ± 0.03 tic sampling of the stellar IMF. This variation is also very small for structures younger than log(τ /yr) = 8.4. Note that the simulated populations younger than this age have typically 10-50 UMS stars. Most observed structures contain a few tens of UMS stars, and some large ones even contain hundreds or thousands of UMS stars. Thus, we conclude that the variation of ξ caused by stochastic sampling is also unimportant.
From the above analysis, we suggest that the stellar number is proportional to the mass for most of the young stellar structures. As a result, the stellar number-size relation can be regarded as the mass-size relation with approximately the same slope. The mass-size relation of substructures in a fractal follows M ∝ R D , where D is, by definition, the fractal dimension (Mandelbrot 1983). Thus, the young stellar structures have a projected fractal dimension of D 2 = κ = 1.48 ± 0.03. We use the subscript "2" to emphasize that the young stellar structures are identified from a two-dimensional KDE map of the projected stellar distributions (Section 4). D 2 reflects significant lumpiness of the young stellar structures, since it is much smaller than 2, the latter being expected for a smooth spatial distribution.
We also investigated the mass-size relations for subsamples of young stellar structures at different significance levels. The corresponding power-law slopes, obtained via least-squares fits, are displayed in Table 2. We find that the mass-size relation has a dependence on the significance levels, since the 1-2σ structures have steeper slopes than those on higher significance levels. This behavior is similar to that found by Gouliermis et al. (2017) for the spiral galaxy NGC 1566. Figure 13 (top panel) shows the size distribution of all identified young stellar structures. The distribution peaks at ∼10 pc, with a drop toward small sizes. This drop is caused at least in part by the incompleteness of young stellar structures (see also Sun et al. 2017a,b). On the one hand, the KDE map (Fig. 7), from which the structures are identified, has a spatial resolution of 10 pc (Section 3). As a result, structures smaller than 10 pc are incomplete due to the finite spatial resolution. On the other hand, we have arbitrarily required each structure to contain at least N min = 5 stars (Section 4). Although this criterion can reject less reliable detections, it may In the bottom panel, the histogram without error bars corresponds to all identified young stellar structures, while the histogram with error bars corresponds to structures larger than R = 10 pc. The two histograms overlap at small Σ, and the latter one is fitted with a lognormal distribution, as shown by the solid-line curve. In all three panels, the error bars reflect Poissonian uncertainties. also remove some genuine structures. This effect can also be assessed with the stellar number-size and surface density-size diagrams (Fig. 11). In the diagrams, the spatial resolution and minimum stellar number are indi-cated by the dashed lines. It is obvious that structures at small R are incomplete, especially for those with low surface densities.

Size, Mass, and Surface Density Distributions
The size distribution becomes increasingly noisy toward large sizes. In the range of 10-100 pc the size distribution is consistent with a single power law. Thus, we have performed such a fit to the data, which is shown as the solid line in Fig. 13 (top panel). The best fitting slope is α(R) = − 1.4 ± 0.1. Substructures inside a fractal follow a cumulative size distribution where D is the fractal dimension (Mandelbrot 1983). Because Eq. 3 is a single power law, it is mathematically equivalent to a differential size distribution (per logarithmic interval, as in Fig. 13) Equation 4 has the same form as the size distribution of young stellar structures with 10 < R < 100 pc. Thus, these structures agree with a fractal dimension of D 2 = −α(R) = 1.4 ± 0.1, which is very close to that derived from the mass-size relation (1.48 ± 0.03; Section 5.3).
At R = 300 pc the power law drops below 1, but there are still a few structures larger than this size. We have checked the boundaries of these structures, whose linear scales are found to be comparable to that of the SMC's main body 8 . As mentioned in Section 4, they trace the extent of the SMC's galactic structure rather than internal structures on sub-galactic scales. We will discuss later (Section 6.1) that the formation of the internal, sub-galactic structures are possibly controlled by supersonic turbulence in the ISM. In contrast, those larger structures are obviously more influenced by global galactic processes. As a result, they do not follow the same power-law size distribution as the sub-galactic ones. Given their small number, we do not attempt to investigate their size distribution in greater detail.
Similar to the size distribution, the stellar number distribution (Fig. 13, middle panel) is also subject to incompleteness at small values. This is evident in the stellar number-size diagram (Fig. 11 top), which clearly shows that the spatial resolution of 10 pc rejects many small-N * structures. This effect is no longer important for structures with N * > 30 stars. Thus, we have fitted the distribution in the range of 30-1000 stars with a single power law. The best-fitting result is shown in Fig. 13 (middle panel), with a slope of α(N * ) = −1.0 ± 0.1. As demonstrated in the previous subsection, N * is proportional to structure mass. Thus, the stellar number distribution reflects the structures' mass distribution, which is also a power law with nearly the same slope. It can be derived mathematically that the mass distribution (per logarithmic interval) of substructures inside a fractal should have a power-law slope of −1, irrespective of the fractal dimension (e.g. Elmegreen & Falgarone 1996, their Section 4). This is in agreement with the measurement obtained here.
The surface density distribution is shown in the bottom panel of Fig. 13. However, it is not easy to assess the influence of incompleteness, since surface density is not correlated with size (Fig. 11, bottom panel) -the effects of spatial resolution and minimum stellar number reject structures of a wide range of surface densities, especially at small sizes and small surface densities. Based on the large scatter of the surface density-size relation, we assume that structures larger than 10 pc should be representative of their underlying young stellar structure population in terms of surface density distributions. Their surface density distribution is shown in Fig. 13 (bottom panel) as the histogram with error bars. The distribution is in good agreement with the solid-line curve, which is a lognormal function fitted to the data.
The appearance of a histogram may depend sensitively on the bin size (Ivezić et al. 2014). We have tested this effect by changing the bin sizes of the parameter distributions. The power-law slopes of the size/stellar number distributions and the lognormal form of the surface density distribution do not show any significant deviations. Thus, the results in this subsection are robust against this effect.
6. DISCUSSION 6.1. Imprints from turbulent ISM There are noticeable similarities between the young stellar structures and the ISM. Similarly to the hierarchical young stellar structures, the ISM also contains significant substructures in a hierarchical manner, including clouds, clumps, and cores on all scales. The ISM substructures exhibit highly irregular morphologies. The projected boundaries of molecular clouds have been investigated based on the perimeter-area relation in a number of studies (e.g., Beech 1987;Scalo 1990;Falgarone et al. 1991;Vogelaar & Wakker 1994;Lee et al. 2016). The typical perimeter-area dimensions are 1.4-1.5, very close to D p = 1.44 ± 0.02 as derived in Section 5.2 for the young stellar structures (although smaller values have also been reported for several molecular clouds; see e.g. Dickman et al. 1990;Hetem & Lepine 1993). Also similar to the young stellar structures, the ISM clouds have power-law size/mass distributions, and they exhibit power-law mass-size relations with fractional slopes (e.g. Elmegreen & Falgarone 1996;Roman-Duval et al. 2010).
Supersonic turbulence has also been argued to play a dominant role in controlling the hierarchical substructures in the ISM (Elmegreen & Scalo 2004). In starforming complexes, supersonic turbulence, aided by gravity, leads to the fragmentation of ISM into smaller and smaller substructures. Stanimirović & Lazarian (2001) analyzed the position-position-velocity data cubes of H I lines in the SMC. They found that the observations agree with theoretical expectations of turbulence (Lazarian & Pogosyan 2000). Thus, the significant density fluctuations of H I in the SMC (Stanimirović et al. 1999(Stanimirović et al. , 2000 are indeed caused by active turbulence instead of being a static structure. With power-spectrum analysis, Stanimirović et al. (1999Stanimirović et al. ( , 2000 found fractal dimensions of D 2 = 1.4-1.5 for H I and dust in the SMC. These values are also consistent with that of the young stellar structures (D 2 = 1.48 ± 0.03 from the mass-size relation, and D 2 = 1.4 ± 0.1 from the size distribution; Sections 5.3 and 5.4).
Thus, the above discussion supports a scenario of hierarchical star formation in a turbulent ISM, which leaves imprints of the ISM properties on the young stellar structures. The remarkable similarities between the ISM and the stellar structures also suggest that the dynamical evolutionary effects are not significant. With dynamical evolution, the young stellar structures in the SMC, especially the unbound ones, will be dispersed on a timescale of ∼100 Myr (see the discussion in Section 3.1).
The effect of turbulence is also consistent with the surface density distribution. Lognormal distributions have been observed for the column/volume densities of molecular clouds (Lombardi et al. 2010;Schneider et al. 2012) or of simulational turbulent gas (Klessen 2000;Federrath et al. 2010;Konstandin et al. 2012). With turbulence, such a distribution can be understood in a purely statistical way, if self-gravity and thermal pressure are unimportant (Vázquez-Semadeni 1994; see also Section 4.4 of Sun et al. 2017b). With star formation following gas distribution, the young stellar structures should also have a lognormal surface density distribution. This is in agreement with the result described in Section 5.4.
Strong self-gravity, shocks, and rarefaction waves may lead to deviations from lognormal density distributions (Klessen 2000;Federrath et al. 2010). Especially, strong gravity in star-forming clouds would cause a power-law tail at high surface densities (see e.g. Schneider et al. 2012). However, our results do not show such deviations to any statistical significance. Such high-surface density structures (e.g. dense star clusters) usually have small sizes. And since our identification is complete only for scales larger than 10 pc, it is possible that the surface density distribution, shown as the histogram with error bars in Fig. 13 (bottom panel), does not reach the high surface densities of the power-law tail. While the all-structure distribution (histogram without error bars) does exhibit a number excess at high surface densities, we note again that this distribution is affected by incompleteness because structures with low surface densities and small sizes could be missing (see also Fig. 11, bottom panel). Gieles et al. (2008) investigated the spatial distribution of stars of various ages in the SMC. Using two-point correlation functions (TPCFs) and the Q-parameter (introduced by Cartwright & Whitworth 2004), they found that the SMC's overall stellar distribution evolves from a high degree of substructures (∼10 Myr) to a smooth radial density profile. At ∼75 Myr the stellar distribution becomes statistically indistinguishable from a smooth distribution, suggesting that most of the young stellar structures have been erased by that age. Bonatto & Bica (2010) studied the size distributions of young (ages 10 Myr) clusters, old (ages 600 Myr) clusters, and non-clusters (nebular complexes and their stellar associations) in the SMC. They were all found to follow power-law size distributions, with slopes of −3.6, −2.5, and −1.9, respectively. Note that their slopes are for the size distributions per linear interval, and correspond to −2.6, −1.5, and −0.9 for the distribu-tions per logarithmic interval. These values differ from α(R) = −1.4 ± 0.1 derived in Section 5.4. On the one hand, star clusters and associations are not distinguished and they are analyzed together in our work. On the other hand, their analysis is based on the catalog of extended sources of Bica et al. (2008), which was compiled from extensive previous catalogs by different authors. As a result, we expect our analysis to be more systematic and less biased. Gouliermis et al. (2014) reported a bimodal clustering of stars in the SMC star-forming region, NGC 346. They found that the TPCF is a double power law with a break at 20 pc. Based on simulations, they argued that the inner, flatter power law reflects a centrally concentrated stellar distribution (typical for star clusters), while the outer, steeper part corresponds to an extended, hierarchical component. Although their study focused on a smaller region in the SMC, they derived a fractal dimension for the hierarchical component of D 2 = 1.4 based on the slope of the outer TPCF (or a volume fractal dimension of D 3 = 2.4 based on their simulations). This value is consistent with what we have derived for the SMC globally.  studied the young stellar structures in the LMC. They found that the structures have a power-law luminosity function (per linear interval) with a slope of 1.94 ± 0.03. Assuming that the structures have no significant age spreads or stochastic sampling effects, this luminosity function corresponds to a structure mass distribution similar to that of the SMC young stellar structures (Section 5.4). They also derived a fractal dimension of D 2 = 1.8 based on the Q-parameter. However, we note that this method is subject to large uncertainties (see Fig. 5 of Cartwright & Whitworth 2004).

Comparison with Other Works
In Sun et al. (2017a,b) we studied young stellar structures in the LMC's 30 Dor and Bar complexes, respectively. Young stellar structures in the two complexes have similar properties to those in the SMC, in terms of the size/mass/surface density distributions and the mass-size relation. More quantitively, the fractal dimensions are D 2 = 1.6 ± 0.3 for the 30 Dor complex and D 2 = 1.5 ± 0.1 for the Bar complex. These values are also close to that for the SMC young stellar structures, irrespective of the different galactic environments for the Magellanic Clouds.
The hierarchical young stellar structures have also been investigated in Galactic star-forming regions and in some other galaxies (e.g. Larson 1995;Simon 1997;Elmegreen & Elmegreen 2001;Elmegreen et al. 2006Elmegreen et al. , 2014Gouliermis et al. 2010Gouliermis et al. , 2015Gouliermis et al. , 2017. The results are in general agreement with ours. Specifically, the typical fractal dimensions are close to 1.4-1.5, but Elmegreen et al. (2014) found that starburst dwarfs or H II galaxies have larger projected fractal dimensions if they contain one or two dominant stellar complexes. We refer to Sun et al. (2017a) for an extensive comparison of the fractal dimensions.

SUMMARY AND CONCLUSIONS
In this paper we study the hierarchical young stellar structures in the SMC. This is based on the VMC survey and UMS stars selected from the (J − K s , K s ) CMD. We apply a contour-based clustering analysis to the KDE map of the UMS stars. Young stellar structures are identified as overdensities on a series of significance levels. After rejecting less reliable detections, we identify 556 young stellar structures in total.
The young stellar structures are distributed in a hierarchical manner, such that larger structures on lower significance levels contain one or more smaller ones on higher significance levels. We illustrate this behavior with a dendrogram, which is a structure tree showing the "parent-child" relations between the young stellar structures on different significance levels. The structures also have highly irregular morphologies. This is quantitively analyzed based on the perimeter-area relation of their (projected) boundaries, which suggests a perimeter-area dimension of D p = 1.44 ± 0.02.
Size, stellar number (the number UMS stars within their boundaries, proportional to mass to a good approximation), and surface density are also analyzed statistically for the young stellar structures. There is a power-law correlation between stellar number and size with a slope of κ = 1.48 ± 0.03, while the surface density and size are not significantly correlated. The young stellar structures follow power-law size and mass distributions (per logarithmic interval) with slopes of α(R) = −1.4 ± 0.1 and α(N * ) = −1.0 ± 0.1, respectively. The surface density of the structures, on the other hand, exhibits a lognormal distribution.
These properties are remarkably similar to those of the ISM which is regulated by supersonic turbulence. Thus, our results support a scenario of hierarchical star formation, which leaves the imprints of the ISM properties on the young stellar structures. We also make a comparison with previous studies of young stellar structures.