THE SPATIAL DISTRIBUTION OF THE YOUNG STELLAR CLUSTERS IN THE STAR-FORMING GALAXY NGC 628

Star formation leads to the creation of stellar clusters (Lada & Lada 2003) and most, if not all, stars form in some type of clustered structure. Observations of local star-forming regions have shown that clustering is a common feature, resulting from the fractal properties of the interstellar medium (ISM) under the effects of turbulence (Elmegreen & Efremov 1997; Elmegreen et al. 2014). As a result, gravitationally bound clusters occupy the smallest and densest regions of the hierarchy of giant molecular cloud (GMC) complexes forming on the larger scales (∼1 kpc).

As clusters can be observed to greater distances than individual stars, young stellar clusters provide an excellent means to investigate the connection between the continuous distribution of star formation on small scales to galactic formation at large scales. While studies on the hierarchical clustering of stellar populations have become more sophisticated over time (e.g., see the early papers of Payne-Gaposchkin 1974; Efremov 1995), whether or not the stars and the clusters are mapping the same type of hierarchy (i.e., the dense peaks are randomly distributed within each hierarchy or they are a biased representation of each hierarchy) is still a question that needs to be answered. The answer will help understand the role of cluster formation and evolution that is complementary to studies of the cluster formation efficiency. The spatial distribution of newly formed stars is also important because it provides insight not only on the processes of star formation, but also on the evolution and environmental dependences within stellar clusters and their migration from their clustered complexes. For instance, we expect that the observed hierarchical clustering will disappear with age (Elmegreen et al. 2006; Elmegreen 2010): the densest regions with the shortest mixing timescales lose their substructures first, whereas the larger, unbound regions will lose their substructure over time owing to random initial motions and tidal forces. The migration timescale for which stars and clusters "abandon" their clustered natal structure is not well constrained for most systems, but such knowledge is vital toward understanding how star formation evolves in both space and time.

One of the most powerful ways to probe the clustering distribution of galactic components is with the two-point correlation function (Peebles 1980; Zhang et al. 2001; Odekon 2006, 2008), a statistical tool used to quantify the excess probability of finding one object within a specified distance of another object against that of a random, unclustered distribution. Applying correlation functions to study the clustering of clusters will provide insight to the (1) physical process of cluster formation; (2) the extent that the formation of clusters is hierarchical; and (3) whether or not the clustered distribution of young stellar clusters reflects the fractal structure of the interstellar gas (Elmegreen & Efremov 1996; Efremov & Elmegreen 1998; Elmegreen & Elmegreen 2001; Bastian et al. 2007).

Thus far, few studies have been done on the clustering distribution of stellar clusters, which shows the need for a systematic approach to the problem across a wide variety of galaxies and environments. In one of the first studies on the hierarchy occurring within stellar cluster ensembles, Zhang et al. (2001) found that a power-law correlation function well described the clustering hierarchy up to 1 kpc for the star clusters within the Antennae galaxies. Scheepmaker et al. (2009), in a study of the stellar clusters in M51, found that clusters showed an age-dependency in their degree of clustering, where the youngest clusters were more clustered compared to older clusters.

Investigations of clustering behavior across time and scale lengths on stellar systems have been performed in a few local galaxies in recent years as well. In a study of the stellar samples of M 31, Gouliermis et al. (2015a) found that the youngest (<25 Myr) stars were more clustered than older stars (<300 Myr) and the observed clustering changed within different regions of the galaxy. Studies of the stellar structures in both the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC) have shown that stars are born highly fractal and evolve toward a uniform distribution within the crossing time of the galaxy (∼175 and 80 Myr for the LMC and SMC, respectively; Gieles et al. 2008; Bastian et al. 2009; Bonatto & Bica 2010). Within NGC 1313, Pellerin et al. (2007) found that infant mortality—the process by which most clusters do not survive the death of the most massive stars in the cluster, forcing the remaining unbound clusters to dissolve over a short period of time ($\lesssim 10\;{\rm{Myr}}:\;$Lada & Lada 2003; Fall et al. 2005; Bastian & Goodwin 2006; Chandar et al. 2010)—is an efficient process in disrupting, and therefore, destroying stellar clusters by showing that B-stars can be found outside star clusters. These unbound stars in the diffuse field can account for a large portion of the UV light from a galaxy (Meurer et al. 1995; Hoopes et al. 2001; Tremonti et al. 2001; Crocker et al. 2015) and the disruption and dissolution of clusters is crucial in populating the diffuse stellar field. The evolution of the clustering of star clusters with time will thus depend on the combination of two factors: how fast the cluster migrates within a galaxy and how long the typical star cluster survives.

If the cluster dispersion timescale has a dependency with galactic environment, we would expect to see their age distribution change over different ambient environments within a galaxy. Indeed, in a study by Sánchez-Gil et al. (2011), the young stellar populations within NGC 628 exhibited age gradients across the spiral arms, with the youngest star formation regions concentrated along the center of the spirals, and an age gradient from the inner to the outer parts of the galaxy. Radial stellar migration can further flatten the observed age gradient, implying steeper age gradients than those measured. Additionally, Bastian et al. (2011a) have found strong dependency between cluster age and galaxy position, where the age distribution varies as a function of galactocentric distance, becoming shallower toward the outskirts of the galaxy. Theoretical work (Elmegreen & Hunter 2010; Kruijssen et al. 2011) suggests that the age distribution of clusters does vary with environment and that it can influence the dissolution timescale, with clusters disrupting faster when the gas surface density ${{\rm{\Sigma }}}_{{\rm{gas}}}$ is high and living longer when ${{\rm{\Sigma }}}_{{\rm{gas}}}$ is low. On the other hand, if cluster dispersion is independent of environment, the age distribution of all clusters should be the same, regardless of position within a galaxy. Studies of both LMC and SMC have shown that there is a flat age distribution of clusters, concluding that cluster disruption does not significantly affect the age distribution, and that there is little dependence on environment for the first ∼200 Myr (Gieles et al. 2007; Baumgardt et al. 2013). Early results of the Panchromatic Hubble Andromeda Treasury (PHAT; Dalcanton et al. 2012) survey show that the clusters of M 31 show no evidence for cluster dissolution at early times (30–100 Myr). Using the first PHAT clusters catalog of 601 clusters, Fouesneau et al. (2014) confirms that cluster disruption has little to no effect prior to the timescale of 100 Myr in Andromeda. The spatial correlation between clusters will allow us to study the extent of evolution in the clusters and whether or not the formation of clusters is hierarchical (Efremov & Elmegreen 1998; Bastian et al. 2007).

Ambient environment is crucial for how clusters form, and quite possibly, it is even more critical to investigate whether stellar clusters trace star formation. Studies on the clustering of stars have revealed important information about the nature of the star formation process itself. The recent results of Gouliermis et al. (2015a) show an environmental dependence in the clustering distribution of stars in the nearby spiral M33, where stars in the outer regions of the galaxy experienced less disruption and larger amounts of clustering compared to stars located in the inner regions. This agrees with the work of Silva-Villa et al. (2014), which showed an environmental dependence on the cluster population within M83, where the clusters in the outer region showed less disruption and flatter age distributions compared to clusters in the inner region. Work by Gouliermis et al. (2015b) for NGC 6503 shows that younger stars follow a hierarchical distribution, whereas older stars display a homogeneous, less clustered distribution. By comparing the clustering results within this work to the clustered results of stars, we will be able to infer if and how star clusters also trace star formation.

The nearby spiral galaxy in this study, NGC 628, was observed as part of the Legacy ExtraGalactic UV Survey²⁰ (LEGUS; Calzetti et al. 2015), a Cycle 21 Hubble Space Telescope (HST) program that has imaged 50 nearby (∼3.5–12 Mpc) galaxies in five UV and optical bands (NUV, U, B, V, and I) with the UV/Visible (UVIS) channel on the Wide Field Camera 3 (WFC3). The aim of LEGUS is to investigate the relation between star formation and its galactic environment in nearby galaxies, over scales ranging from individual star systems to kpc-sized structures. NGC 628 is a face-on grand design spiral galaxy with a large number of star clusters to provide a statistically powerful test for changes in clustering strength with scale. The occurrence of hierarchical structures within NGC 628 was already investigated by Elmegreen et al. (2006) using ACS data in B, V, I, and Hα bands. These authors examined the distributions of size and luminosity of star-forming regions, finding that both were well described with a power law, indicative of a hierarchical structure of stellar components within the galaxy.

The work herein primarily focuses on describing the two-point correlation function as a tool to study the clustering properties of the young stellar clusters. For this goal, we use the galaxy NGC 628 as an example. We hope to address the following points once the clustering properties of a larger number of LEGUS galaxies have been studied: (1) whether star clusters are clustered and how strongly; (2) how the clustered distribution of stars compares to that of star clusters; (3) how the clustering depends on age and environment and the dissolution timescale of the clusters; and (4) whether or not clusters can be used to trace the structure of star formation.

The galaxy selection and cluster identification process is described in Section 2. The methodology of the two-point correlation function is introduced in Section 4. In Section 5, we describe the results and analysis and how we use the correlation function to draw conclusions about the properties of our star clusters. We discuss our results concerning the hierarchy of the stellar clusters in Section 6. Finally, we summarize the findings of this study in Section 7.

We study the face-on grand design spiral galaxy NGC 628 located at a distance of 9.9 Mpc (Olivares et al. 2010) with no apparent bulge. NGC 628 is the largest galaxy of the M74 galaxy group and has a global SFR(UV) of approximately 3.7 ${M}_{\odot }\;{{\rm{yr}}}^{-1}$ (Lee et al. 2009). The stellar disk appears largely undisturbed in studies of the gas kinematics with optical and UV imaging (Herbert-Fort et al. 2010) and has a disk thickness of 0.25 kpc (Peng 1988). We select this galaxy primarily because of the low-inclination angle (25°), relatively large angular size (105 × 95), clearly defined spiral arms, and high number of clusters available for analysis.

The images of NGC 628 for our analysis were obtained in five broadband filters: NUV (F275W), U (F336W), B (F438W), V (F555W), and I (F814W) in 2013 October with the UV-optical channel (UVIS). Observations of NGC 628 consist of two pointings: one toward the center of the galaxy, NGC 628c, and one toward the galaxy's eastern region, NGC 628e, adjacent to the first pointing. All five filters in both pointings are aligned and referenced to the UVIS/F336W filter. Both pointings are combined into a single mosaic for analysis, shown in Figure 1. The images have a scale of 0.04 arcsec per pixel, corresponding to a pixel resolution scale of 1.9 pc at a distance of 9.9 Mpc. General descriptions of the standard data reduction of the LEGUS data sets is available in Calzetti et al. (2015).

Zoom In Zoom Out Reset image size

Stellar clusters within NGC 628 are identified by first undergoing an automated process using SExtractor (Bertin & Arnouts 1996) that avoids background sources and single, bright stars. Each catalog includes sources that satisfy the following two conditions: (1) the V-band concentration index (CI) must be greater than the stellar CI peak value (CI = 1.4 and 1.3 for the inner and outer pointing, respectively, in order to make the selection criteria more uniform between the varying pixel scales of the cameras); and (2) the source should be detected in two contiguous bands with a signal-to-noise greater than three.

As stars appear unresolved even at the highest HST resolution power, their CI will vary little and their CI distribution will be highly picked around an average value typical of a stellar point-spread function. On the other hand, clusters are partially resolved and their sizes can vary, therefore, they generally have a CI larger than that of stellar values. Using the CI distributions of the extracted cluster candidates, we selected the smallest CI value that would allow us to remove the bulk of the stellar interlopers from our catalog. Because the resolution power is also dependent of the HST camera used, the values used for the inner and outer frame of NGC 628 are different. Simulations of CI as function of effective radius of the clusters show us that the chosen limits correspond to ${R}_{{\rm{eff}}}=1$ pc at a distance of NGC 628.

In order to secure reliable measurements from spectral energy distribution fitting, each source is required to have a 3σ detection in at least four of the five photometric bands, which is necessary to adequately break the age-extinction degeneracy. The physical properties of each cluster (ages, extinctions, masses) have been derived using deterministic stellar population models (Yggdrasil; Zackrisson et al. 2011) and an ${\chi }^{2}$ fitting approach as described in Calzetti et al. (2015), which includes uncertainty estimates (see Adamo et al. 2010). In this paper we use standard catalogs that contain cluster physical properties derived with deterministic models with solar metallicity for both stars and gas, an average covering factor of 50%, and a Milky Way extinction curve (Cardelli et al. 1989) with a foreground E(B – V) = 0.06. Photometry is performed with a circular aperture of 4 pixels in radius, with the background measured within an annulus of 7 pixels in inner radius and 1 pixel in width. Single stellar populations are used to determine the ages of our clusters; for clusters that are better described with multiple stellar populations, our reported age is recovering the mean age of the cluster.

The uncertainties we derive in age, mass, and extinction are about 0.1 dex. Within the LEGUS collaboration we also undertake cluster analysis based on stochastically sampled stellar libraries and Bayesian approaches (Krumholz et al. 2015). A comparison between the standard approach based on Yggdrasil deterministic models (Zackrisson et al. 2011) and the analysis performed with SLUG models (da Silva et al. 2012) shows that the average derived cluster properties do agree down to cluster masses of $5000\;{M}_{\odot }.$ Below this mass range we observe the largest uncertainties in derived ages, masses, and extinctions because, as widely discussed in the literature (e.g., Cerviño & Luridiana 2004), deterministic models are not able to correctly interpret the large scatter in color caused by the stochastic sampling of the cluster IMF. However, in this analysis we are mainly focusing on a morphological division of the cluster catalogs based on visual inspection. We also explore the change in clustering using two age ranges (i.e., 1–40 and larger than 40 Myr) and two mass ranges (below and above log(M) = 3.6 ${M}_{\odot }$). A realistic error to associate to the value of 40 Myr is the average uncertainty of 0.1 dex observed in deterministic cluster properties; this corresponds to 40(+20/−10) Myr. The low mass range will be affected by stochastic effects and we discuss this point in Section 4.3. More details will be given in a forthcoming paper (A. Adamo et al. 2015, in preparation).

Each object within the cluster catalog that has an absolute magnitude brighter than −6 mag (well within the 90% completeness limits) in the V band is then visually inspected and assigned one of four categories: (1) a symmetric, centrally concentrated cluster; (2) a concentrated cluster with some degree of asymmetry; (3) a multiple peaked system; or (4) a spurious detection, such as a foreground/background source, single bright star, bad pixel, or a source that lies too close to the edge of the chip. Class 4 objects are excluded from the final cluster catalog. Figure 2 shows what a typical cluster looks like for each classification. Classifications 1, 2, and 3 are considered to be genuine star clusters or associations. The classification of each object is compiled using visual inspection from at least four independent members of the LEGUS team. The final cluster catalog is compiled by comparing all the results from each individual, after checking for consistency.

Zoom In Zoom Out Reset image size

The distribution of masses and ages of the clusters is shown in Figure 3, with the three classes identified by different colors. In the same figure, the selection limit of ${M}_{{\rm{F}}555{\rm{W}}}=-6.0$ mag (apparent magnitude of 23.98 mag) for the clusters is shown by a continuous line, which follows the characteristic age-mass correlation for ages $\gt {10}^{7}\;{\rm{year}}.$ General properties of the three cluster classes that will be detailed below include that Class 3 clusters tend to be on average younger and less massive than Class 1 and Class 2 clusters.

Zoom In Zoom Out Reset image size

Cluster candidates that meet the requirements of a minimum cutoff of CI $\gt \;1.4$ (CI $\gt \;1.3$ for the clusters in NGC 628e) and detection in at least four filters but are fainter than an absolute V-band magnitude of −6 are not visually classified; however, these unclassified sources are still added to the cluster catalogs, labeled as Class 0 cluster candidates. We are greatly interested in using the faint cluster candidates to increase our sample, especially at the low mass end (Krumholz et al. 2015). Because these are not visually inspected sources, we acknowledge that there will be contamination of spurious objects in this class. We impose additional constraints on the Class 0 (not visually inspected) objects by analyzing those with CI $\leqslant$ 1.6 mag separately from those with CI > 1.6 mag. The additional criterion of selecting the broadest of the Class 0 sources will increase our chance of identifying genuine clusters within the classification (see Section 5.2 for a more detailed discussion). For all Class 0 sources, we place a limit on the broadest CI of 4 to exclude background galaxies that may be included in the catalog. We identified 345 Class 0 cluster candidates. For the 128 Class 0 candidates with CI > 1.6 mag, we performed a visual check to quantify the accuracy of these objects being clusters. The results can be found in the appendix; we find that the contamination rate of non-clusters in Class 0 is as high as 57%.

NGC 628c(e) has a total of 1019 (245) star clusters classified as Class 1, 2, or 3 and 345 (39) Class 0 cluster candidates with more than four filters available. We combine both pointings for our analysis, which results in a total of 1264 genuine (Class 1, 2, and 3) and 384 Class 0 cluster candidates. There are 17 clusters in the overlapping region of the east and central pointings that are present in both catalogs; we removed the duplicate clusters within the NGC 628e catalog. The breakdown of all the classifications in age, mass, CI, and color excess E(B – V) are shown in Figure 4. Figure 5 shows the breakdown of only Class 0 objects, with a CI cut at 1.6 mag, where we expect the sources with the broadest CI indices to exclude all non-cluster objects. As expected, the clusters that are classified as multiple peak sources (Class 3) have the largest CI value on average. Class 1 (centrally concentrated) clusters have larger masses and older ages when compared to the rest of the clusters in the galaxy. The stellar clusters (Class 1, 2, 3) show small extinction, with an average color excess of E(B – V) ∼ 0.15. On the other hand, Class 0 candidates exhibit moderate extinction, with an average of E(B – V) = 0.46, which remains roughly the same for all Class 0 clusters, independent of CI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Correlation functions provide statistical means by which to measure and provide detailed quantification on the clustering distributions of galaxy constituents. One of the most commonly used quantitative measures of clustering structure is the three-dimensional two-point correlation function, $\xi (r),$ which measures the magnitude of clustering as a function of scale size.

Originally defined by Peebles (1980) as ${dP}={\bar{n}}^{2}[1+\xi (r)]{{dV}}_{1}{{dV}}_{2}$ for cosmological applications, the two-point correlation function $\xi (r)$ is defined as a measure of the probability of finding a neighboring object, above what is expected for an unclustered random Poisson distribution, in a shell element with volume element ${dV}$ at a distance r from any object with mean number density of the galactic structure, $\bar{n}.$ The two-point correlation function, $\xi (r),$ is usually fit as a power law, $\xi {(r)=(r/{r}_{0})}^{-\gamma },$ where r₀ is the characteristic scale length of the clustering, defined as the scale at which $\xi (r)=1.$ The correlation function has been observed to have the same slope γ across galaxy and cluster systems (e.g., Peebles 1975; Postman et al. 1998; Daddi et al. 2000; Brodwin et al. 2008); the uniformity of the measured indices seems to suggest a common underlying dynamic on all scales.

The spatial distribution of galactic components can also be measured in two dimensions as projected onto the plane of the sky. In our study, we measure the two-dimensional (2D) projected angular correlation function $\omega (\theta ),$ defined as the probability above Poisson of finding two star clusters with an angular separation θ as ${dP}={N}^{2}[1+\omega (\theta )]\ d{{\rm{\Omega }}}_{1}d{{\rm{\Omega }}}_{2},$ where N is the surface density of clusters per steradian with two infinitesimal elements of solid angle $d{{\rm{\Omega }}}_{1}$ and $d{{\rm{\Omega }}}_{2},$ separated by angle θ.

To measure $\omega (\theta ),$ pairs of stellar clusters are counted as a function of separation, compared by what is expected for an unclustered distribution. A random catalog (in x, y position) of sources must be produced, populating the same sky coverage and geometry (e.g., edges, masks) as the data with randomly distributed points. We define masks as areas that exclude all data, such as the ACS chip gap, or a reduction in the observed surface area of available data with respect to the global average, such as dust lanes or interarm regions of the galaxy. We account for masks in the random catalog that are present in the real data. The ratio of pairs of clusters observed in the data relative to pairs of points in the random catalog is then used to estimate $\omega (\theta ).$ In this study, we implement the Landy & Szalay (1993, LS) estimator, which cross-correlates the data and random catalog to minimize edge and mask effects, written as

$$\begin{eqnarray}&&{\omega }_{{\rm{LS}}}(\theta )=\displaystyle \frac{\mathrm{DD}(\theta )-2\mathrm{DR}(\theta )+\mathrm{RR}(\theta )}{\mathrm{RR}(\theta )},\end{eqnarray} \tag{ 1 }$$

where DD is the number of data–data pairs, DR is the number of cross-correlated data–random pairs, and RR is the number of random–random pairs with the same mean density and sampling geometry with separation between θ and $\theta +\delta \theta .$ The pairs are computed as

$$\begin{eqnarray}\mathrm{DD}(\theta ) & = & \displaystyle \frac{{P}_{\mathrm{DD}}(\theta )}{N(N-1)}\\ \mathrm{DR}(\theta ) & = & \displaystyle \frac{{P}_{\mathrm{DR}}(\theta )}{{{NN}}_{R}}\\ \mathrm{RR}(\theta ) & = & \displaystyle \frac{{P}_{\mathrm{RR}}(\theta )}{{N}_{R}({N}_{R}-1)},\end{eqnarray} \tag{ 2 }$$

where N and N_R are the total number of data and random points in the survey volume, respectively, and ${P}_{\mathrm{DD}}(\theta ),\;$ ${P}_{\mathrm{DR}}(\theta ),$ and ${P}_{\mathrm{RR}}(\theta )$ represent the total pair counts in each separation $\theta \pm \delta \theta$ bin for the data–data, data–random, and random–random pairs, respectively. The size of the bin is determined by the sample size; we calculate our correlation function with eight spatial bins, spaced logarithmically between 016 (4 pixels) and 200'' (5000 pixels), corresponding to spatial scales of 4.8 pc to 9.6 kpc, which is the largest scale we can investigate with our data. Bin numbers were selected as a compromise between resolution and the total number of clusters available to be sampled within each bin. When computed as above, a random distribution of an unclustered population will result in a flat correlation $\omega (\theta )=0,$ while a clustered distribution will have $\omega (\theta )\gt 0.$ For hierarchical structuring, a general trend of decreasing $1+\omega (\theta )$ with radius is expected (Gomez et al. 1993; Larson 1995). This correlation function has also been successfully used for characterizing the clustering behavior of both stars and stellar clusters (e.g., Gomez et al. 1993; Zhang et al. 2001; Scheepmaker et al. 2009; Gouliermis et al. 2015a, 2015b) in several nearby galaxies. Although the formulation in Equation (1) is optimized for taking into account edge effects, we still attempt to reproduce as closely as possible the geometry of the galaxy region sampled. Because NGC 628 is nearly face-on, we do not need to take into account deprojecting the data before computing the correlation function.

It is also important that the random catalog must be large enough to not introduce Poisson error in the estimator. This is checked by ensuring that the RR pair counts in the smallest bin are high enough, such that Poisson errors are negligible in the total error budget. As a result, the number of random points needs to be greater than the data points. We achieve this for a random catalog that is 100 times greater than the size of the cluster sample. The uncertainties in the measurements of $\omega (\theta )$ in each radius bin are estimated as a Poisson error.

4.1. Application of the Clustering Estimator

Applying the two-point correlation function to our star clusters helps us identify common age structures and derive the correlation length as a function of age and location. Figure 6 shows the correlation function for both pointings of NGC 628 as well as the mosaic. It is interesting to note that the Class 1 (centrally concentrated) clusters have a nearly flat relationship across all spatial distributions, while the strongest clustering is seen in Class 3 (multiple peak) sources. The correlation for Class 0, which potentially may also include spurious sources, is flatter than that for our most structured cluster but still has a distribution that is more fractal (i.e., steeper exponent) than Class 1. In Sections 5.1 and 5.2 we discuss the differences between the correlations for the separate classes of clusters, and how we attempt to single out the Class 0 clusters without visual inspection.

Zoom In Zoom Out Reset image size

As mentioned above, the LS estimator is minimally sensitive to the random sample size and is effective at handling both edge and masking effects, as found in Landy & Szalay (1993) and Kerscher et al. (2000). Using the central pointing, NGC 628c, to test how changing the geometry and masks within the random data sample is going to affect the results, we found that a random catalog of the same observing geometry as the real data with a random sample size of 10,000 objects already converges to the same correlation function values as a random catalog with 500,000 sources. As a sanity check, we create random fields with 1200 and 70,000 objects and compute the correlation function between the two catalogs and find that the computed $\omega (\theta )$ is zero across all angular scales.

We also verify that the LS correlation function is only minimally sensitive to the presence of masks within the data. After populating the random catalog for NGC 628c, we first removed the small chip gap (see Figure 7) of several hundred pixels that runs across the length of the detector from the random catalog. In order to test how the missing data alters the correlation function, we simulated dust lanes in the random catalog that coincide with the spiral arms in the galaxy where there is a dearth of observed star clusters. Adding anywhere between a single dust lane (about 6% of the total area) to up to four dust lanes (13% of the total data area) has a negligible effect on the resulting correlation function, shown in Figure 8. However, the clusters are already undersampled in the interarm regions and their removal along with the random objects in the same region, as expected, does not affect the results. Removing the random sources from these undersampled regions results in an expected decrease of the observed clustering at small scales because it removes the data–random cross-correlation count pairs for the clusters that reside next to these regions, but does not have other major effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To test the effect of masks in regions that are not undersampled, we take a random catalog with the same geometry as covered by the real data, remove regions from the random catalog that cover 10% or 20% of the total area within the galaxy, and recompute the correlation function. We keep all the data and just remove the region from the random field. We perform this test 10 times, moving the region of missing data around the galaxy, and take the average resulting correlation function, as shown in Figure 9. For areas that are outside the dust lanes (i.e., we are not selecting regions where the clusters are particularly undersampled), we observe a slight decrease in the clustering at small spatial scales and a transition toward increased clustering at larger spatial scales. Because we are removing regions that are populated more than the dust lane regions, this can be best understood as an extreme example of removing data in regions that are not undersampled. The increase of clustering strength at the largest spatial scales is similar to what is observed when we artificially increase the random catalog size. These differences are not particularly significant from the statistical point of view.

Zoom In Zoom Out Reset image size

We finally test the case of a random catalog with area that is 10% or 20% larger than that of the actual data (Figure 10). For this test, we normalize the number of objects within the random catalog with the size of the random field geometry so that each field has the same surface density of random sources. In this case, we observe an increase in the strength of the clustering (higher values for the two-point correlation function). Likewise, a random catalog that has a smaller area than the coverage of the real data will result in a decrease in clustering across all scales. The effect of matching area and shape becomes especially important for galaxies that do not span the entire width of the chip and have to be described with ellipses for a subset of the chip coverage. While not applicable to NGC 628, it does affect other LEGUS galaxies in the survey.

Zoom In Zoom Out Reset image size

5.1. Quantifying the Correlation Strength

The projected angular two-point correlation function (Equation (1)) is usually well described with a power law:

$$\begin{eqnarray}&&\omega (\theta )={A}_{\omega }{\theta }^{\alpha },\end{eqnarray} \tag{ 3 }$$

where the slope α measures the strength of the clustering and the amplitude A measures the correlation length of the clustering; we use both to determine if the clustering is consistent with being scale-free. In a fully hierarchical model, the clusters will be correlated with other clusters from other groupings on larger scales, together forming a group higher up in the hierarchy of star formation. This would result in a smooth decline of $1+\omega (\theta )$ in both radius and age. Figure 11 shows the two-point correlation function after combining both pointings in NGC 628; we show each cluster class in a separate panel, and fit each correlation function with a broken power law (except for Class 0 and 1, which are better fit with a single power law). We compute the fit and determine the slopes using the Levenberg–Marquardt nonlinear least-square minimization fit, where the functional form is given as

$$\begin{eqnarray}&&\mathrm{log}[1+\omega (\theta )]\\ \quad =\left\{\begin{array}{lr}{A}_{1}+{\alpha }_{1}\mathrm{log}(\theta ) & :\mathrm{log}(\theta )\lt \beta \\ {A}_{2}+({\alpha }_{1}-{\alpha }_{2})\beta +{\alpha }_{2}\mathrm{log}(\theta ) & :\mathrm{log}(\theta )\gt \beta ,\end{array}\right.\end{eqnarray} \tag{ 4 }$$

and where the break point β is the logarithm of the position of the separation break along the x-axis (spatial scale), A₁ and A₂ are the clustering amplitudes before and after the break, and ${\alpha }_{1}$ and ${\alpha }_{2}$ are the slopes to the power law before and after the break point, respectively. Slopes, amplitudes, and break-points are free parameters in the fit. The fitted parameter results for NGC 628 are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Power-law Parameters of NGC 628

Class	Number	A1	${\alpha }_{1}$	β	A2	${\alpha }_{2}$
				(arcsec)
Class 1	420	2.24(3)	−0.14(3)	...	...	...
Class 2	432	5.31(3)	−0.82(10)	3.3	2.76(2)	−0.17(2)
Class 3	412	8.17(7)	−1.51(13)	1.2	3.4(2)	−0.37(5)
Class 0	384	3.44(5)	−0.20(2)	...	...	...
Class 0; CI > 1.6	128	4.8(1.2)	−0.31(5)	...	...	...
Class 1, 2, 3	1264	3.91(19)	−0.65(9)	3.3	2.29(2)	−0.18(2)
Class 1, 2, 3; CI > 1.6	575	3.9(3)	−0.75(5)	3.3	2.26(2)	−0.161(14)

Note. Columns list the (1) Classification of stellar clusters—Class 1: symmetrical clusters, Class 2: asymmetrical clusters, Class 3: multiple peak clusters, and Class 0: non-visually identified clusters; (2) Number in each classification; (3) Amplitude A₁ of the angular correlation function before the break point; (4) Slope ${\alpha }_{1}$ of the angular correlation function after the break point; (5) Location of the break point β—cluster classifications that are best fit with a single power law do not have a break point; (6) Amplitude A₂ of the angular correlation function past the break point; and (7) Slope ${\alpha }_{2}$ of the angular correlation function after the break point. Numbers in parentheses indicate uncertainties in the final digit(s) of listed quantities, when available.

Download table as:  ASCIITypeset image

For a hierarchical (scale-free or self-similar) distribution, the total number of clusters N within an aperture increase with radius r as $N\propto {r}^{\alpha +2},$ where α is related to the (two-dimensional or projected) correlation fractal dimension D2 $=\;\alpha +2$ (e.g., see Falgarone et al. 1991). The distribution of star formation and interstellar gas over a large range of environments is observed to show a (three-dimensional, 3D) fractal dimension of D3 ∼ 2.3 (e.g., Elmegreen & Efremov 1996; Elmegreen & Elmegreen 2001), in agreement with the predicted fractal dimension for the density structure of a turbulent ISM (Federrath et al. 2009). Projection of this 3D fractal on a plane will result in a 2D fractal dimension D2 = D3 − 1 if the perimeter–area dimension of a projected 3D structure is the same as the perimeter–area dimension of a slice of the structure (Elmegreen & Scalo 2004). The two-point correlation function results are qualitatively identical to a power spectra analysis. However, as shown in simulations by Gouliermis et al. (2014), the D2/D3 conversion is not trivial.

As can be seen in the power-law fits to NGC 628 in Figure 11, the power-law index recovered ranges from $\alpha =-0.14$ to −1.51, covering spatial scales from 5 pc to nearly 5 kpc. If we take the weighted mean of the measured indices before the break point at 33, we recover an average slope of $\alpha \sim -0.8,$ which is in agreement with what is expected from a hierarchical distribution of gas. After the break point, we observe a dramatic decrease in the fractal distribution, informing us that the spatial distribution of the clusters (excluding Class 1 and 0 sources) is systematically less clustered at larger spatial length-scales. Results of a study of the hierarchical structure of star formation within a subset of 12 LEGUS galaxies by Elmegreen et al. (2014), suggests that hierarchically structured star-forming regions are common unit structures, which can be several hundred parsecs, but also that the observed self-similarity of young stellar structures down to parsec scales is indicative that individual star clusters form stellar groupings on larger scales.

Given the 27 field of view of WFC3/UVIS, the largest physical size that can be probed for a single pointing in NGC 628 is 7.8 kpc. A pixel scale corresponds to a physical size of 1.90 pc at the distance of 9.9 Mpc, with a typical star cluster around 5–10 pc in diameter. Each of the classes is best fit with a broken power law at 33, corresponding to a spatial scale length of 158 pc at a distance of 9.9 Mpc, except for Class 1 clusters, best fit with a single power law. The break in the power law corresponds to the observed line of sight thickness of the galactic disk (see Sections 5.1.2 and 6 for discussion). Class 3 clusters exhibit the break below the rest of the clusters at 12 in addition to exhibiting a much steeper slope across all scale lengths. The exact location of the break point naturally depends on the number of bins used to calculate the correlation function (See Section 3). Past the 33 break, all the classes show a correlation function with a power-law index of $\alpha \sim -0.2,$ indicating that any occurring clustering has been considerably weakened by this length scale, lasting until a few kiloparsecs. The clustering observed increases below this break point, albeit to varying degrees between the cluster classifications. These results indicate that star formation is clumpy and the stellar clusters form highly clustered distributions for spatial separations of r less than the break point. For spatial scales above the break point, the clustering observed becomes more homogeneous and less clustered.

In order to address how clusters behave versus associations, we compare Class 1 + 2 (bound) clusters to just Class 3 clusters in Figure 12. We believe that Class 1 and 2 are gravitationally bound clusters in different states of relaxation, whereas Class 3 (multiple peak) clusters are quite possibly  unbound associations. We see that the clustering behavior of the Class 1 + 2 is different than what is observed for the Class 3 clusters alone. These results indicate that bound clusters (Class 1 and 2) have a correlation behavior that is weaker, corresponding to a different spatial distribution when compared to what we observe for Class 3 and the compilation of Class 1, 2, and 3 clusters combined. This strengthens our results that as (bound) clusters relax, they become less clustered and behave differently than what is observed for associations.

Zoom In Zoom Out Reset image size

5.1.1. Break in the Power Law

Excluding Class 0 and 1 sources, distributions are described with a broken power law with the break point occurring at 158 pc averaged between all the genuine clusters (Class 1, 2, and 3), although Class 3 clusters lose their self-correlation at much shorter scale lengths—58 pc. In fact, we may expect a smaller characteristic scale for the younger clusters and a larger one for older clusters, if clusters and stars represent the same hierarchy (Efremov 1995). Class 1 clusters are the least self-correlated of all the cluster morphologies with a nearly flat exponent, implying that the distribution of these clusters is nearly uniform across all scales. The flattening of the exponent found after the break point for Classes 2 and 3 is comparable to the uniform distribution of Class 1, implying the loss of self-correlation distribution for these types of clusters after the break point at the larger scale. For the cumulative correlation function of Class 1, 2, and 3, we still see strong self-correlation with a break point below 158 pc, suggesting that the Class 2 and 3 clusters are dominating the small-scale correlation.

A similar break in power spectra within other galaxies has been observed and interpreted as due to the line of sight thickness of galactic disks (Lazarian & Pogosyan 2000; Elmegreen et al. 2001; Odekon 2008; Block et al. 2010) when the line of sight depth is smaller than the transverse size. A study of the brightest stars in the LMC, M31, and M33 by Odekon (2008) identified a transition in the correlation function to a higher correlation dimension (weaker clustering), with the transition marking the large-scale regime where disk geometry and dynamics set the scale for structure. This transition from the small-scale 3D turbulence to large-scale two-dimensional disturbances in the disk, which are observable with a transition to a more shallow power law at the transition, has been detected for the HI gas in the LMC, with a transition occurring at 180 pc by Padoan et al. (2001) and Elmegreen et al. (2001), and 290 pc by Kim & Park (2007). We may be seeing the analogous transition within our young stellar populations.

5.1.2. Dependency of Results on Bin Numbers

The choice of annuli width used to calculate the correlation function in Equation (1) will naturally influence the results and the exact location of the break point of the correlation behavior. In order to ensure adequate sampling of clusters in each bin across the scale length that we are investigating, our bin numbers were selected such that the smallest angular bin contained a minimum of two clusters, because the smallest separation size is where our data are the sparsest. Figure 13 shows how changing the bin size from five bins to 20 bins moves the location of the break point, averaged over Class 1, 2, and 3 clusters, between 2 and 14 arsec (96 pc and 672 pc, respectively). Decreasing the bin numbers aids to smooth out occurring variations at both small and large angular separation. We can see that the break point only starts to disappear at the coarsest sampling of 5 bins, highlighting that galaxies with a large number of clusters available are needed to provide statistically sound results. In order to see the break point at the disk thickness, we need to be able to resolve the equivalent scale.

Zoom In Zoom Out Reset image size

5.2. Age and CI

The correlation function of Class 0 objects tends to be very shallow, confirming that Class 0 objects contain a compilation of every possible object available. In order to attempt to identify only the star clusters within Class 0 objects, we take very conservative cuts in CI and then further divide into age bins. The youngest objects that are genuine star clusters should exhibit an increase the strength of the power law in the correlation function. We see this in the correlation function for Class 3 clusters, which on average are younger and more clustered compared to Class 1 clusters. In order to test if the correlation changes with cluster age, we divide our clusters by age, starting at 10 Myr, and recompute the correlation function. Figure 14 shows the correlation function for all cluster classifications as well as only the clusters within each classification that are younger than 40 Myr, which is the age where we see the most noticeable transition from clustered to non-clustered. Excluding Class 0 and 1 clusters, the correlation function only increases when we consider the youngest clusters in each class. Separating the broadest sources (CI > 1.6 mag) for the combination of the Class 1, 2, and 3 sources from their narrower counterparts (far right, bottom panel of Figure 14) only marginally increases the clustering strength. The shallowness of the slope for the Class 1 sources is not entirely surprising because these are the oldest and most massive clusters in the galaxy, where we may be viewing a subset of clusters that have dispersed and lost their clustered, nonhomogeneous distribution. Class 3 sources, which are on average the youngest in the galaxy, are systematically the most clustered class observed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen in our analysis, at 20 Myr we see a transition at which the clustering structure disappears and the clusters display a more homogeneous, non-clustered distribution, in place by ages of 40 Myr. However, the clustering strength of Class 3 clusters is only minimally effected with an increase in age. All star clusters regardless of their classification tend to be randomly distributed above this age, where random motions and shear effects can explain the observed data trends. The flatness of the correlation function of the age distribution for Class 1 clusters as compared to Class 3 clusters can be evidence of cluster dissolution because we see a dramatic decrease in the total strength of the clustering with increasing age.

We also show the Class 0 sources in Figure 14 divided into CI bins above and below the CI of 1.6 mag, where the broader sources (CI > 1.6 mag) do show evidence of increased clustering. We do not show the Class 0 sources for ages older than 40 Myr because there is a dearth of available data. When we combine the youngest and broadest sources, limiting our selection to our very conservative cut of CI $\geqslant 1.6$ mag and divided into extremely young ($\leqslant 10$ Myr) and young ($\leqslant 40$ Myr) age bins, we do not see a difference in the strength of the clustering. However, out of all the CI $\geqslant 1.6$ mag sources (128), 98% have ages of $\leqslant 40\;{\rm{Myr}},$ which indicates that sources with high values of CIs are sources with younger ages.

As we saw in the power-law fits of Section 5.1, Class 0 sources are described with a fairly shallow power law of index $\alpha =-0.20,$ which is not as steep as expected if the objects are part of a larger hierarchical distribution. However, the clustering does increase when we limit our selection to Class 0 clusters with CI > 1.6, as can be seen in Figure 14. If we repeat the power-law fits to the CI $\gt \;1.6$ clusters, we find that the power-law index increases to $\alpha =-0.31,$ up from $\alpha =-0.20.$ This is still not as steep as the slope observed for Class 2 or 3 clusters, however, it compares to slopes measured by other star cluster studies in the range of $-0.4\ {\rm{to}}\;-0.8$ (Scheepmaker et al. 2009) and $-0.7\ {\rm{to}}-1$ (Zhang et al. 2001). Figure 15 shows that while the correlation steepens for the broader (CI > 1.6) sources in Class O, however, the age has no influence on the strength of the correlation function.

The age where we start to see a decrease in the clustering of stellar clusters, 40 Myr, is significantly younger than the age at which the clustering structure of stars starts to dissipate (300 Myr, e.g., Bastian et al. 2011b). The difference in age can be explained if clusters lose their substructure through a combination of both dissolution and migration. While stars maintain their individuality over their lifetime, clusters can dissolve or evaporate over a much shorter timescale, implying that a cluster may lose its identity while the stars that compose it do not.

5.3. Mass

We investigate the role that the mass of stellar clusters plays in the hierarchy, formation, and dispersion of stellar structures.

Figure 16 shows how the clustering changes when we divide the clusters in the sample into bins of high and low mass; the value $\mathrm{log}(M/{M}_{\odot })=3.6$ was chosen to ensure adequate numbers of clusters sampling each bin across all cluster classifications. For Classes 1 and 2, we notice slightly stronger correlations for the lower mass clusters, where the correlation strength starts to increase for the lower mass clusters at the smallest spatial scales. The strength of the clustering for different masses is not significantly different for Class 3 clusters. Combining all the genuine clusters (bottom middle and right plots in Figure 16) accentuates the increased clustering strength for lower mass clusters. We do recognize that clusters with $\mathrm{log}\lt 3.6\ {M}_{\odot }$ may be biased by stochastic IMF sampling, resulting in unreliable mass measurements at the low mass end. This makes it difficult to draw specific conclusions on how the mass affects clustering and is further compounded by the fact that we lack clusters populating the lowest spatial scales on the correlation function. For the Class 0 sources that are limited to CI > 1.6, a large increase in clustering is observed for the large mass bin of $\mathrm{log}(M/{M}_{\odot })\gt 3.6.$ If we consider all the Class 0 sources regardless of CI value, this same jump is not seen and there is no significant difference between the low and high mass sources. This increase in clustering for the more massive clusters is not seen for any of the other cluster types within NGC 628.

Zoom In Zoom Out Reset image size

The hierarchical morphology in stellar ensembles has been investigated by Elmegreen et al. (2014) with 12 LEGUS galaxies, including NGC 628. The power laws of star-forming regions suggest that hierarchically structured regions are represented with a common unit length of a few hundred parsecs. The observed power law of the hierarchical structure of star-forming complexes observed in galaxies is consistent with the model that star formation is regulated by turbulence, where gas compression forms successively smaller clouds within larger clouds (see Vázquez-Semadeni et al. 2009; Elmegreen 2010; Kritsuk et al. 2013; Elmegreen et al. 2014). The hierarchical nature of star formation will likewise drive the observed hierarchy of young stars, and therefore, the distributions of young stellar clusters. The same processes will likely form a secondary correlation for age, making larger regions older in proportion to the turbulent crossing time (Efremov & Elmegreen 1998). These results indicate that the young stellar clusters we observe should also be imprinted with the hierarchical structure of the natal cloud structure from which they are born. The influence of the self-similar nature of star formation is demonstrated on the distributions of resolved massive young stars on galactic scales and the power laws in their parameter correlations, as shown for NGC 6822 (Gouliermis et al. 2010), as well as with the much deeper LEGUS data for NGC 6503 (Gouliermis et al. 2015b). Our results are consistent with the study of the luminosity function of the young stellar populations in NGC 628 by A. Adamo et al. (2015, in preparation), which were found to be well described by a power-law distribution with index close to −2, suggesting that clusters form in a turbulence-driven fractal ISM.

Excluding Class 0 and 1 sources, our correlation functions for NGC 628 are best described with two power laws connected at a break point corresponding to 158 pc, although this break is much shorter for Class 3 clusters. The thickness of the galactic disk of NGC 628 is around 0.25 kpc (Peng 1988; Ma et al. 1998), a factor of two different from our measured break point in the power spectrum. Previous results of power spectra in galactic observations have shown that breaks in the power law are related to the line of sight thickness of galactic disks. Observations of the LMC disk with IR emission by Block et al. (2010) show a two-component power spectrum with a break at 100–200 pc, occurring at a depth that is comparable to the disk line of sight thickness. The same break in the power spectra for the LMC at 100 pc is also seen with HI observations (Elmegreen et al. 2001), interpreted as the line of sight thickness. If the line of sight thickness is responsible for the observed break in power spectra, dwarf galaxies could lack the break as the line of sight depth is comparable to the transverse length (Westerlund 1997; Roychowdhury et al. 2010). Increasing the scale height results in an increase of chance alignment along the line of sight, serving to decrease both the clustering scale length and magnitude of clustering.

Interestingly enough, the youngest population ($\leqslant 10$ Myr) of star clusters observed by Scheepmaker et al. (2009) was observed to have the shallowest slope ($\alpha =-0.4$), where their oldest star clusters ($30\leqslant {\rm{Age}}\;{\rm{(Myr)}}\lt 400$) have the steepest slope of $\alpha =-0.8.$ Our youngest clusters on average (Class 3) display the most prominent clustering behavior, whereas the oldest clusters on average (Class 1) systematically show minimal clustering with a very shallow slope. Despite being fairly flat, the clustering results of our young clusters within Class 0, which have not been visually classified and could still include non-clusters, appear to be consistent with other observations, even if they fall short of the measured slopes of $\alpha \sim -0.8$ by Zhang et al. (2001). While these early results seem to suggest that the correlation function of young clusters may be similar between different galaxies (i.e., it hints toward a universal fractal dimension of hierarchical star formation), more studies of different galaxies are necessary to determine the validity of this statement.

As shown in Figure 14, the strength of the clustering decreases with both increasing spatial scale and increasing age, excluding our multiple peak clusters (Class 3), although the increase in the clustering is only marginal for the symmetric clusters (Class 1). At ages younger than 20 Myr, we see a transition from the stellar clusters displaying a highly clustered distribution toward a flat, non-clustered distribution for ages greater than 40 Myr. This is consistent with star clusters born within a highly clustered structure and the clustering dissipating in a little as a few Myr after the formation from random motions of the clusters and shear effects.

We see a slight increase in clustering with decreasing mass (Figure 16), which is predominatly only affecting the smallest spatial scales. This is not an expected result as in a hierarchy/turbulence model universe; everything is self-similar so there should be no mass dependency on the clustering results. We know that there is an age effect on the clustering correlation, and to understand whether there is a bias present that makes our selection of low mass clusters younger, we plot clusters younger than 40 Myr divided into high and low mass bins, as shown in Figure 17. When only considering our youngest clusters, the dependency of the mass of a cluster on the clustering strength disappears and we can see that the correlation is nearly identical for clusters above and below our mass bin of 3900 M_⊙. Thus the mass dependence seen in Figure 16 is probably an age effect (i.e., low mass clusters tend to be younger and more correlated than high mass clusters). We lack number statistics in NGC 628 to be able to investigate how mass influences the clustering behavior for clusters older than 40 Myr and we defer it to future studies.

Zoom In Zoom Out Reset image size

We see that the opposite behavior occurs for the broadest Class 0 candidate cluster sources because there is a steep rise in the clustering for massive clusters (Figures 16 and 17). However, since we are dealing with small number statistics (there are only two objects within the first two bins of the CI > 1.6 and mass >3900 M_⊙ bin) the clustering strength increase is accompanied with large error bars, making any conclusion unreliable at best. With greater cluster numbers, we will be able to statistically investigate whether the properties of the star clusters change as a function of environment and if this is reflected in a change of the observed clustering properties.

In this paper, we present a technique to investigate the spatial clustering of the young stellar clusters in the LEGUS galaxy NGC 628 with UV+optical data taken with WFC3/UVIS from HST. The inclusion of NUV observations within this study provides reliable measurements of both the age and masses of the clusters, giving us an unprecedented, high-angular resolution view of the clustering in a way that has not been possible before.

Through visual inspection, we identified 1264 stellar clusters within the galaxy that have an absolute V-band magnitude brighter than −6, CI $\gt 1.4,$ and detection in at least four filters. For the fainter (V-band magnitude fainter than −6) clusters that did not undergo visual inspection (a total of 384 objects) but still have detection in at least four filters, we take cuts in both CI and age to separate out possible cluster candidates from spurious, non-cluster candidates. Using this method, we identify 128 possible clusters of out a total of 384 faint sources, increasing the number of possible stellar clusters available for analysis to 1392.

We implement the two-point correlation function to describe and study the amount of clustering occurring within our cluster candidates. We find that the clusters of NGC 628 have correlation functions that can be well described with a power-law slope consistent with what is expected if the clustering is part of a hierarchy and that the observed clustering strength increases within younger clusters. Our results reveal the internal hierarchical morphology present in the spatial distribution of the stellar clusters. We find that the general observed hierarchy of the young stellar clusters decreases monotonically across the dynamical range of the galactic disk from 7 to 10 kpc.

Our clusters with centrally concentrated light profiles have the highest average mass and oldest ages, and show the flattest power-law slope of $\alpha \sim -0.14$ across all spatial scales. Our multiple peaked clusters have the steepest slope of $\alpha \sim -1.51$ and correspond to the clusters with the lowest average masses and youngest ages in the galaxy. The asymmetrical clusters display a slope of $\alpha \sim -0.82.$ The weighted average of the measured power-law slope for the combination of all the genuine clusters is $\alpha \sim -0.8.$ The results from the correlation function reveal that our cluster population is part of a self-similar distribution with an average 2D fractal distribution D2  =  1.2. This fractal distribution is consistent with the hierarchical morphology of star-forming structures that was formed from fractal gas that has undergone turbulent fragmentation.

Both the asymmetrical and multiple peaked sources show a break in their power law, at 33 and 12 respectively, corresponding to a spatial scale length of 158 pc and 58 pc; the average slope at length-scales beyond the break point is $\alpha \sim -0.24.$ These results suggest that the stellar clusters do form in clustered structures at the small spatial scales, while the clustering is erased and the spatial distribution is observed to be more homogeneous at larger spatial scales beyond the break at 33. Conversely, the symmetrical clusters have a single shallow slope of $\alpha \sim -0.14,$ indicating nearly uniform distribution across all scales. The measured slope for the faint non-visually classified sources has a power-law slope of $\alpha \sim -0.20;$ when we exclude sources with measured CI values less than 1.6, the slope steepens to $\alpha \sim -0.31,$ which is still below what we find for the average of the genuine clusters in the survey. The asymmetrical and multiple peak clusters dominate the observed clustering structure.

In addition to the clustering strength decreasing with increasing spatial scale, we also find that the clustering strength decreases with increasing age, which is consistent with the fractal distribution being erased as the clusters age. We find a dramatic decrease in the clustering at a timescale of 40 Myr. We also find that dividing up clusters by mass has a minor influence on the clustering behavior, however, the low mass bin ($\lt 4000$ $\;{M}_{\odot }$) is subject to stochastic IMF effects and as such, these masses are highly uncertain. The mass dependency on the clustering behavior vanishes when we only consider the young clusters below 40 Myr.

Future investigations of a larger sample of LEGUS galaxies will provide a better understanding of whether there is a common correlation length observed in the stellar clusters across galaxies and if the clustering depends on the ambient galactic environment. The clustering nature of the young stellar clusters will not only provide a window to investigate whether or not star formation occurs in hierarchical patterns in both space and time, it will also inform us about the clustering properties of stars at the time of their formation, providing important information on the nature of the star formation process itself. If stars and clusters map the same type of hierarchy, this would allow the use of observations on the clustering distribution of stellar clusters to trace and study star formation, providing a way to understand the cluster formation efficiency in relation to the star formation efficiency.

We thank the anonymous referee for the careful and constructive comments that greatly improved this paper.

This paper is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, under NASA Contract NAS 5-26555. These observations are associated with Program 13364 (LEGUS). Support for Program 13364 was provided by NASA through a grant from the Space Telescope Science Institute.

This research has made use of the NASA/IPAC Extragalactic Database (NED) operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

D.A.G. acknowledges financial support by the German Research Foundation (DFG) through grant GO 1659/3-2. M.F. acknowledges support by the Science and Technology Facilities Council (grant number ST/L00075X/1).

In order to ensure that our method selecting the broad (CI > 1.6) sources within the Class 0 candidates accurately identifies genuine star clusters, we visually checked the 128 Class 0 sources with CI > 1.6, finding that 73 are true star clusters (57%). While this is a high contamination rate and does affect our cluster results, placing a limit for sources at CI > 1.6 does improve the blind selection of clusters and allows us to identify a correlation at a much higher degree of confidence than without the selection (dotted orange versus black in Figure 18). Additionally, small numbers limit our results of NGC 628 and we will only be able to fully investigate more significant results when we employ the galaxies within a larger selection of the LEGUS sample.

Zoom In Zoom Out Reset image size