collapse. Young stars and star clusters inherit the spatial properties of the natal gas from which they form, and can therefore be used as tracers to understand the physical mechanisms at play in the star formation cycle. Unlike young stars though, young stellar clusters can be observed to greater distances, and hence provide an excellent source of information to investigate the complex mechanisms of star formation in diverse environments.

Star formation in galaxies is spatially structured in a hierarchical, scale-free pattern such that, smaller and denser associations that extend all the way down to substellar scales, are surrounded by larger, less dense ones, that go out to kiloparsec scales (see, Elmegreen 2010, for a review). These scale-free structures are analogous to geometric fractals (Mandelbrot 1982), and have been shown to be present in the distribution of unbound stars (see, Gouliermis 2018, and references therein), embedded stars in clusters and starforming regions (Sánchez et al. 2007;Fernandes, Gregorio-Hetem & Hetem 2012;Gregorio-Hetem et al. 2015;Sun et al. 2017), H II regions (Feitzinger & Galinski 1987;Sánchez & Alfaro 2008), OB associations (Bresolin et al. 1998;Pietrzyński et al. 2001;Kumar, Kamath & Davis 2004;Gutermuth et al. 2008), star-forming regions (Elmegreen & Elmegreen 2001;Elmegreen et al. 2006Elmegreen et al. , 2014;;Bastian et al. 2007;Rodríguez, Baume & Feinstein 2020;Mondal et al. 2021), and young star clusters (Zhang, Fall & Whitmore 2001;Bastian et al. 2005;Scheepmaker et al. 2009;Grasha et al. 2015Grasha et al. , 2017a,b),b), and are expected to originate from the inherently hierarchically structured interstellar gas distribution (Elmegreen & Salzer 1999;Elmegreen et al. 2003;Elmegreen & Scalo 2004;Elmegreen 2007;Bergin & Tafalla 2007;Dutta et al. 2009;Federrath, Klessen & Schmidt 2009;Beattie et al. 2019). This hierarchical nature of gas is consistent with that set by the scale-free physical mechanisms that act on it, i.e. gravity and interstellar turbulence, and is central to the so-called gravoturbulent fragmentation (Elmegreen 1993;Klessen, Heitsch & Mac Low 2000;Mac Low & Klessen 2004;Padoan et al. 2014;Federrath 2018) and global hierarchical collapse (Vázquez-Semadeni et al. 2009;Vázquez-Semadeni, González-Samaniego & Colín 2017) scenarios, both leading theories describing the multi-scale star and cluster formation process in the interstellar medium (ISM; McKee & Ostriker 2007;Krause et al. 2020). Individual stars form at the smallest scales of this hierarchy and group together to form star clusters, which themselves are spatially correlated with other star clusters in kpc-scale star complexes and flocculent spiral arms (Efremov 1995;Elmegreen & Efremov 1996;Gusev 2002;Bastian et al. 2005;Ivanov 2005). The spatial correlations of star clusters would thus trace the largest scales of this hierarchy that has properties which would presumably be set by the physical mechanisms governing star formation at galactic scales. Studying the spatial structure of star clusters is thus an effective way of obtaining insights into the physical mechanisms at play.

The two-point correlation function (TPCF) is a robust tool to probe the scale dependence of the clustering properties of a distribution, as it quantifies how much excess correlation a distribution of points has at a given separation (angular or linear) compared to a completely random distribution (Peebles 1980). For a hierarchical or scale-free distribution, a general trend of the TPCF decreasing with separation is expected (Gomez et al. 1993;Larson 1995;Bate, Clarke & McCaughrean 1998). Such a trend has been seen in earlier studies of star clusters in galaxies such as the Antennae (Zhang et al. 2001), M51 (Bastian et al. 2005;Scheepmaker et al. 2009), and NGC 0628 (Grasha et al. 2015). Apart from providing an estimate for the strength of the correlation at a given scale, the TPCF can also quantify the spatial heterogeneity of the distribution through an effective fractal dimension, which can be obtained from the slope of the TPCF (see e.g. Calzetti, Giavalisco & Ruffini 1989;Falgarone, Phillips & Walker 1991). The fractal dimension quantifies how space-filling or clumpy a distribution is, and is expected to be set by turbulence in the ISM (Stutzki et al. 1998;Elmegreen & Scalo 2004;Sánchez, Alfaro & Pérez 2005;Federrath et al. 2009).

It has often been argued that the fractal dimension observed in the ISM has a nearly universal value of around ∼2.3 (Elmegreen & Falgarone 1996), suggesting a universal nature of the self-similar hierarchy. However, more recent work has questioned this univer-sality, especially at galactic spatial scales, finding variations in the inferred fractal dimensions and the scales to which the hierarchy extends. Such differences could arise due to the different sources of turbulence that might dominate at various scales and environments, and set different density structures (see; Federrath et al. 2009), and/or from the modification of the scale-free behaviour due to galactic scale dynamical processes such as rotation, shear, or feedback (Padoan et al. 2001;Odekon 2008;Sánchez et al. 2010;Dib et al. 2020). For instance, Sánchez & Alfaro (2008) investigated the fractal nature of H II regions in 93 nearby galaxies, and found statistically significant variations among the galaxies of the sample, with signs of higher fractal dimensions for brighter, more massive galaxies. Similarly, Grasha et al. (2017a), through the use of the angular TPCF of star clusters in six galaxies observed as part of Legacy ExtraGalactic UV Survey (LEGUS), found a large range of fractal dimensions and correlation lengths, with hints at systematic variations with the galaxy stellar mass and star formation rate. If these results are confirmed, it would suggest that the environment of the host galaxy is important in setting the hierarchical structure of star formation at galactic scales, which would qualitatively be consistent with recent evidence for the same in observations (see recent reviews by Adamo 2015;Chevance et al. 2020) and numerical simulations (Kruijssen et al. 2011;Renaud 2018;Pfeffer et al. 2019).

The spatial distribution of star clusters is also expected to evolve with age. For instance, there is strong evidence that hierarchical clustering dissipates with age, which manifests as a reduction in the TPCF for older populations of stars and stellar clusters (Bastian et al. 2005;Gieles, Bastian & Ercolano 2008;Scheepmaker et al. 2009;Sánchez & Alfaro 2009, 2010;Gouliermis, Hony & Klessen 2014;Gouliermis et al. 2015b;Grasha et al. 2015Grasha et al. , 2017a)). This decrease in spatial correlation implies that older stars/clusters are more randomly positioned than younger ones. In addition, clusters closer to each other tend to have about the same age, regardless of that age, and this leads to an age-difference versus separation relation in the distribution (see e.g. Elmegreen & Efremov 1996;Efremov & Elmegreen 1998; de la Fuente Marcos & de la Fuente Marcos 2009; Grasha et al. 2017b). The randomization of the cluster distribution could be the result of ballistic motion of mutually unbound clusters away from their birth sites, a product of larger scale effects such as shear or tidal interactions, or due to the superposition of successive generations of star formation (see, for instance, Elmegreen 2018, for a discussion). Regardless of the mechanism, the observed time-scale over which the distribution randomizes is ∼40-100 Myr (Grasha et al. 2017a). However, this time-scale also varies from galaxy to galaxy (Grasha et al. 2018(Grasha et al. , 2019) ) and the environment within a galaxy (Silva-Villa et al. 2014;Gouliermis et al. 2015a), and is at least qualitatively consistent with theoretical and numerical work (Elmegreen & Hunter 2010;Kruijssen et al. 2011;Reina-Campos & Kruijssen 2017).

In this study, we use the angular TPCF to investigate the environmental dependencies and evolutionary effects of the hierarchical distribution of star clusters in 12 nearby galaxies as part of the LEGUS (Calzetti et al. 2015a). LEGUS is a Cycle 21 Hubble Space Telescope (HST) Treasury program that imaged 50 nearby (∼3-18 Mpc) galaxies in UV and optical bands, and used the data to identify and prepare catalogs of individual star clusters in the galaxies. We compute and compare the TPCF of the cataloged star clusters among 12 galaxies drawn from the LEGUS sample and search for correlations between the TPCF and the physical conditions of the host galaxy. In addition, we use the estimated ages of the star clusters to probe the time evolution of hierarchical structuring. This study extends the work by Grasha et al. (2017a) to six more galaxies (4): Inclination angle in degrees. References for adopted inclinations in order of the rows: Lang et al. (2020), Koribalski et al. (2018), Lang et al. (2020), Meidt et al. (2009), Lang et al. (2020), Oh et al. (2015), Hunter et al. (1998), Colombo et al. (2014), Koribalski et al. (2018), Walter et al. (2008), Greisen, Spekkens & van Moorsel (2009), Koribalski et al. (2018). ( 5): Position angle measured anticlockwise from the celestial north. References for adopted angles identical to the inclination angles, except for NGC 3738 which adopts the value reported in Vaduvescu et al. (2005). ( 6): Redshift-independent distances adopted from Anand et al. (2021) for all galaxies except NGC 3344, NGC 3738, NGC 5253, and NGC 6503 for which we use the value reported in Sabbi et al. (2018). ( 7): Galaxy integrated Far-UV calculated star formation rate adopted from Calzetti et al. (2015a). ( 8): Stellar mass adopted from Calzetti et al. (2015a). ( 9): Standard isophotal radius of the galaxy adopted from de Vaucouleurs et al. (1991) after applying the distances reported in Column 6. ( 10): Star formation rate surface density obtained by -(i) averaging SFR UV uniformly in a disc of radius R 25 if the entire radial extent of the star-forming gas is contained in the LEGUS field of view, (ii) averaging the local dust-extinction corrected SFR UV values only in the LEGUS field of view, if not. Values in this case are obtained from (Adamo et al., in preparation). ( 11): HI gas surface density obtained by averaging the total HI mass in the disc M HI uniformly in a disc of radius R 25 . M HI values are adopted from Calzetti et al. (2015a). ( 12): Total number of identified star clusters in the LEGUS catalog that we use.

(12 in total), yielding a larger sample size with more statistical power to constrain any potential dependencies on galaxy properties. In addition, we develop a set of physically motivated toy models to interpret the various features we find in our TPCFs, similar to the approach in Gouliermis et al. (2014). The paper is outlined as follows: In Section 2, we describe the galaxies in our study, and, briefly, the procedure adopted by LEGUS to prepare the star cluster catalogs for them. The methodology we adopt to compute the TPCF is provided in Section 3. We present our observed TPCFs and their evolutionary changes in Section 4 along with details on the statistical tools and toy models we adopt to understand the features in the TPCF. Following this, in Section 4.4, we calculate some key physical quantities that highlight differences in the overall spatial distribution of clusters among the galaxies, and attempt to compare this with the properties of the host galaxy. Finally, we summarize our findings in Section 5.

In this study, we select twelve local (<18 Mpc) galaxies from the LEGUS survey of various morphological types, ranging from irregular dwarfs to grand design spirals. The twelve galaxies were picked from the larger sample of LEGUS galaxies for which cluster catalogs were available based on the conditions that -(i) they contain sufficient number of star clusters to calculate TPCFs across a range of separations, and (ii) were relatively face-on to prevent line-of-sight inclination effects. The galaxies and their average physical properties are listed in Table 1, and we provide more detail on the individual galaxies below.

The LEGUS sample consists of both archival and new imaging with either the Wide Field Camera 3 (WFC3) or the Advanced Camera for Surveys (ACS). The F275W (UV) and F336W (U) filters of each galaxy in this sample are WFC3 imaging. The three other bands are taken with either the ACS (archival) or WFC3 (new observations for the LEGUS Programme GO-13364 Calzetti et al. 2015a): ACS/WFC3 F435W (B), ACS/WFC3 F555W/F606W (V), and ACS/WFC3 F814W (I). In this study, we refer to the passbands by the conventional Johnson passband naming: UV, U, B, V, and I, where the V-band is adopted as the reference frame. LEGUS photometry is in the Vega magnitude system. The frames are aligned and rotated with North up. Some of the galaxies are observed with more than one pointing and combined into a single mosaic, whereas others are observed with a single pointing. Reduced science frames in all filters have been drizzled to a common scale resolution, corresponding to the native WFC3 pixel size (0. 03962/px).

NGC 0628 is a nearly face-on spiral galaxy (morphology SAc) located at a distance of ∼9.8 Mpc (Anand et al. 2021). It was observed by the LEGUS survey with two pointings. This galaxy has been observed extensively by all recent major surveys of interstellar gas and dust in nearby galaxies, including THINGS, HERACLES, SINGS, KINGFISH, EMPIRE, and with ALMA (Kennicutt et al. 2003(Kennicutt et al. , 2011;;Walter et al. 2008;Leroy et al. 2009;Bigiel et al. 2016;Turner et al. 2019). Elmegreen et al. (2006) investigated hierarchical star formation in this galaxy, and Grasha et al. (2015) report a measurement of its TPCF.

NGC 1313 is a mildly inclined barred galaxy (morphology SBd) that may be interacting with a satellite, producing a loop of HI gas around the galaxy (Peters et al. 1994) and a recent increase of the SFR in the south-west arm (Silva-Villa & Larsen 2012). Due to both its physical and morphological properties, including the presence of a bar and an irregular appearance, NGC 1313 has been compared to the Large Magellanic Cloud (de Vaucouleurs 1963). Hannon et al. (2019) and Messa et al. (2021) have recently analysed the properties of NGC 1313's star clusters and H II regions, using H α narrow-band and near-infrared Paβ observations. LEGUS observed NGC 1313 with two distinct pointings, and we use a mosaic prepared from the two pointings for the analysis in this study.

NGC 1566, the brightest member of the Dorado group, is an almost face-on spiral galaxy with an intermediate-strength bar and open, knotty arms, a small bulge, and an outer pseudo-ring made from arms that wind antiparallel to the bar ends (Buta et al. 2015). Salo et al. (2010) propose that the spiral arms are formed through bar-driven spiral density waves, and Shabani et al. (2018) find evidence for an age gradient in the star clusters across the spiral arms consistent with the stationary density wave theory. Grasha et al. (2017a) and Gouliermis et al. (2017) use LEGUS catalogs to study the hierarchical distribution of star clusters and young stellar populations, respectively. We caution that there is significant uncertainty in the distance to the galaxy, with published estimates varying from 5.5 to 21.3 Mpc (Tully 1988;Mathewson, Ford & Buchhorn 1992;Willick et al. 1997;Theureau et al. 2007;Tully et al. 2013;Sorce et al. 2014;Calzetti et al. 2015a;Sabbi et al. 2018;Anand et al. 2021). Here, we adopt the value obtained from the Kourkchi-Tully group catalog (Kourkchi & Tully 2017), i.e. 17.7 Mpc, which uses a distance to the Dorado galaxy group obtained through numerical modeling of its orbits. That said, the only effect of changing the adopted distance for our study is that it would shift the correlation functions we obtain along the linear distance axis, and the corresponding conversion from angular separation to linear separation. This galaxy is observed with a single pointing by LEGUS.

NGC 3344 is an isolated barred spiral galaxy (morphology SABbc) with two ring-like morphological features at 1 and 7 kpc, and a small bar within the inner ring (see, for e.g. Verdes-Montenegro, Bosma & Athanassoula 2000). This galaxy was included in the sample studied in Grasha et al. (2017a). Meidt, Rand & Merrifield (2009) analysed the spiral structure and dynamics in this galaxy using H I and CO gas. The LEGUS survey observed NGC 3344 with a single pointing.

NGC 3627 is a strongly barred spiral galaxy (Buta et al. 2015) that has been studied in large atomic and molecular gas surveys (Walter et al. 2008;Leroy et al. 2009;Kennicutt et al. 2011) and exhibits strong burst signatures at the two interfaces of bar and arm in the north and south (Kennicutt et al. 2011), making it a prime candidate for the study of bar-arm interactions (Beuther et al. 2017). In addition, NGC 3627 appears to be interacting with the neighbouring galaxy NGC 3628 (see e.g. Soida et al. 2001), which is expected to be the cause of a perturbed morphology of its western arm, and a higher H 2 /H I mass ratio relative to other local star-forming galaxies (e.g. Saintonge et al. 2011). LEGUS observed this galaxy in a single pointing.

NGC 3738 is an irregular dwarf galaxy (morphology Im) in the Messier 81 group classified as a blue compact dwarf (BCD). It is close to the Milky Way (∼9.9 Mpc), and has a relatively small size (R 25 ∼ 4 kpc). Hunter et al. (2012) included this galaxy in the LITTLE THINGS HI survey, and reported that the HI component of NGC 3738 is morphologically and kinematically disturbed, possibly a result of an advanced merger or ram pressure stripping (Ashley et al. 2017). In addition, Hunter et al. (2018) found that the gas pressure, density, and star formation rate in a localized region in the south-west part of the galaxy is much higher than the rest of the galaxy, with a higher fraction of younger clusters found there. LEGUS observes the entire extent of the optical galaxy in a single pointing.

NGC 4449 is an irregular barred starburst galaxy (morphology SBm) with ongoing and intense star formation distributed along a barlike structure with two streams stemming from its ends. The gas component shows morphological features that may be caused by dynamical interactions with neighbouring galaxies (Hunter et al. 1998). It has a rich population of young, intermediate, and old star clusters, making it the best sampled dwarf galaxy in our study, potentially due to a rich star formation history sculpted by earlier interactions and mergers (see for e.g. Cignoni et al. 2019). Whitmore et al. (2020) include this galaxy in the recent H α-LEGUS survey that add narrowband H α imaging to a subsample of LEGUS galaxies, allowing the production of new cluster catalogs with improved ages. However, since age accuracy is not a significant constraint for our analysis, we choose to use the original LEGUS catalogs for consistency with the remainder of the sample.

NGC 5194 (M51a or the Whirlpool galaxy) is a well-studied spiral galaxy (morphology SAbc) due to its large size, relative proximity, and almost face-on inclination. This galaxy contains the largest number of clusters in our sample. Its grand-design morphology, high-star formation rate, rich-star formation history, and numerous star-forming complexes and star clusters make it a benchmark for nearby extragalactic surveys (e.g. PAWS Schinnerer et al. 2013). A number of authors have investigated NGC 5194's TPCF (Bastian et al. 2005;Scheepmaker et al. 2009) as well as cross-correlations between clusters and molecular clouds (Grasha et al. 2019). In addition, Messa et al. (2018a,b), as part of the LEGUS survey, study the age and mass distributions of the young star cluster population and their dependencies on the local environment within the galaxy. NGC 5194 is known to be interacting with its companion galaxy (NGC 5195) resulting in a marked spiral geometry along with a tidal tail due to the interaction, and the two galaxies together are referred to as the M51 system. We note that the LEGUS field of view was obtained through multiple pointings, and the cataloged star clusters cover both members of the system. However, we found that removing the contribution from NGC 5195 star clusters does not change the TPCF, since it contains a very small fraction of the overall catalog, and hence we keep the overall catalog for completeness. We also refer to the system simply as NGC 5194 since this is the major contributor to the observed TPCF, and doing so maintains consistency in the format of names for the galaxies in this study.

NGC 5253 is a nearby BCD galaxy that hosts a very young central starburst, likely triggered by infalling material along the minor axis of the galaxy (Meier, Turner & Beck 2002;Turner et al. 2015;Miura et al. 2015Miura et al. , 2018)). This results in a dense, clumpy, central region, hosting a rich population of dense super star clusters with very high-star formation efficiencies (Turner & Beck 2004;Calzetti et al. 2015b;Turner et al. 2017;Smith et al. 2020). This galaxy has the lowest number of cataloged clusters in this study, with an age distribution that skews young. The radial extent of the star cluster population is entirely covered in the LEGUS field of view, observed with a single pointing.

NGC 5457, commonly referred to as the Pinwheel Galaxy, is a relatively large, almost face-on (i ∼ 18 • ) SABcd-type spiral galaxy with a complicated arm structure and a highly asymmetric disc morphology suggestive of previous accretion or interaction (Waller et al. 1997;van der Hulst & Sancisi 1988;Walter et al. 2008). It has 823 cataloged star clusters spread across the extent of its large disc. Due to its large angular size, LEGUS covers this galaxy with five different pointings: one in the central region, three in the north-west, and one in the south-east. In this study, we only use the star clusters in the region spanned by the available mosaic galaxy image available on the LEGUS public website, 1 since our analysis method requires knowledge of the mosaic footprint (see below). This mosaic does not span the south-east regions of the galaxy, and hence we exclude these star clusters from our TPCF analysis. Overall, this results in a relatively lower completeness in the azimuthal and radial sampling of the young clusters as compared to the smaller galaxies in our sample.

NGC 6503 is a spiral galaxy classified as an SAcd type in de Vaucouleurs et al. (1991). It has well-developed spiral arms and traces of a bar (Buta et al. 2015). The galaxy has a patchy circum-nuclear appearance in the gas distribution, a morphology that carries over to the young star cluster distribution observed with LEGUS. Freeland et al. (2010) interpreted this as an inner ring around the galactic bar. Gouliermis et al. (2015b) studied the hierarchical distribution of unbound stars with the LEGUS survey, and found that younger stars are organized in a distribution with a 2D fractal dimension of 1.7, whereas older stars display a homogeneous distribution, with a structure dispersion time-scale of ∼ 60 Myr. NGC 6503 also has a significant line-of-sight inclination, for which we compensate by de-projecting star cluster positions before computing the TPCF. This galaxy is fully covered with a single pointing.

1 https://archive.stsci.edu/prepds/legus/dataproducts-public.html 2.1.12 NGC 7793 NGC 7793 is a flocculent spiral galaxy (morphology SAd). It is part of the Sculptor group, and is one of the closest galaxies in the LEGUS sample. It is characterized by diffuse, broken spiral arms with no bar, a very faint central bulge, and a relatively low star formation rate. Sacchi et al. (2019) study the star formation history of this galaxy using LEGUS data, while Grasha et al. (2017a,b) study the spatial and temporal TPCF of its star clusters. In addition, Grasha et al. (2018) study the connection between molecular clouds and young star clusters, and the time-scales of their mutual association in this galaxy. LEGUS observed this galaxy with two pointings, one each in the eastern and western parts of the galaxy.

A detailed description of the standard data reduction of the LEGUS sample can be found in Calzetti et al. (2015a) and in-depth descriptions of the cluster extraction, classification, photometry, and SED fitting procedure are detailed in Adamo et al. (2017). The procedure to obtain the catalogs for the dwarf galaxies are given in Cook et al. (2019). We refer the reader to these papers and provide a brief description of the LEGUS cluster catalogs here.

Catalog construction in LEGUS is a multistep process that begins with an initial automated extraction of cluster candidates identified using Source Extractor (SEXTRACTOR; Bertin & Arnouts 1996) from the white-light images produced with the five standard LEGUS bands (Calzetti et al. 2015a). The SEXTRACTOR parameters are optimized to extract sources with at least a 3σ detection in a minimum of five contiguous pixels.

The automatic catalogs for each galaxy includes sources that satisfy the two following conditions: (1) the V-band concentration index (CI ≡ magnitude difference of a source in an aperture of 1 pixel compared to an aperture of 3 pixels) must be greater than the stellar CI peak value; and (2) the source must be detected in at least two contiguous filters (the reference V band and either B or I band) with a photometric error σ λ ≤ 0.35. These conditions minimize stellar contamination and yield cluster candidates with a signal to noise greater than 3, which allows for reliable constraints on the derived cluster properties of age and mass. This procedure produces our automated cluster catalog that is complete for clusters down to 1 pc in size for galaxies at distances up to 10 Mpc (Adamo et al. 2017). This size is well below the peak of the size distribution of star clusters of ∼3 pc (Ryon et al. 2017).

The next step in catalog construction is photometry. The analysis pipeline measures the luminosity of each cluster using a science aperture of radius 4-6 pixels depending on the distance to the galaxy, with sky corrections computed using a sky annulus at 7 pixels with a width of 1 pixel. The pipeline then applies an average aperture correction, which it estimates from the difference between the luminosity within the science aperture and that within a 20-pixel aperture with a 1-pixel sky annulus for a control sample of isolated clusters. The pipeline applies this aperture correction independently in each filter. The photometry reported in the catalog is in the Vega magnitude system and is also corrected for foreground Galactic extinction (Schlafly & Finkbeiner 2011).

To ensure that we can derive reliable estimates of cluster physical properties (age, mass, and extinction), the next analysis step is to remove from the catalog any clusters that lack a 3σ detection in at least four of the five photometric bands. The pipeline then estimates the masses and ages of the remaining clusters by fitting the observed SED using Yggdrasil deterministic stellar population models (Zackrisson et al. 2011) using a χfoot_1 fitting approach that includes uncertainty estimates (Adamo et al. 2010(Adamo et al. , 2012)). The uncertainties derived in the physical parameters for the final LEGUS star clusters are on average 0.1 dex (Adamo et al. 2017). The Yggdrasil models are based on STARBURST99 (Leitherer et al. 1999) stellar population spectra coupled with nebular emission computed using CLOUDY (Ferland et al. 1998(Ferland et al. , 2013)). Adamo et al. (2010) adopt a Kroupa (2001) IMF in the range 0.1-120 M (see, however, Ashworth et al. 2017, for a generalization to a variable IMF). We also adopt the Adamo et al. (2017) catalogs that use the Padova stellar isochrones that include thermally pulsating asymptotic giant branch stars (Girardi et al. 2000;Vázquez & Leitherer 2005) and a starburst attenuation curve (Calzetti et al. 2000) with the assumption that stars and gas undergo the same amount of reddening. Adamo et al. (2017) computes the nebular emission assuming a hydrogen number density n H = 10 2 cm -2 , a covering factor c = 0.5 (i.e. 50 per cent of the Lyman continuum photons that are produced by the central source driving nebular emission from within the LEGUS aperture), and a gas filling factor f of 0.01, typical of H II regions (Croxall et al. 2016).

The final step in the analysis is visual classification. Members of the LEGUS team visually examine cluster candidates in the automated cluster catalog if they satisfy the following two criteria: (1) detection in a minimum of four bands (VBI and U and/or UV) with a S/N above 3 sigma and (2) brighter than -6 mag in the V-band (Grasha et al. 2015;Adamo et al. 2017). The cluster catalog for NGC 5194 is obtained with a combination of visual and Machine Learning procedures for the final cluster classifications (Grasha et al. 2019). The human or machine classifiers assign each cluster to one of four morphological classes: Class 1 contains compact, symmetric, and centrally concentrated clusters. Class 2 includes compact clusters with asymmetry. Class 3 are compact associations that show multiple-peaked profiles on top of an underlying diffuse emission. Class 4 is the label given to non-cluster contaminants that remain in the catalog after all selection criteria. These are usually bad pixels, foreground stars, or background galaxies. Star cluster candidates that are not visually inspected are labeled as Class 0 in the final catalogs. In this study, we consider all cluster candidates (class 1, 2, and 3) for our results and analysis and do not separate by cluster classification type. We note that some of the candidates we chose to represent as 'star clusters', particularly class 3 ones, may not be gravitationally bound and could disperse or dissolve in relatively short time-scales (∼10 Myr). However, since we are interested in the hierarchy of star formation and the spatial distribution of star formation, including these unbound associations in our definition of star clusters is warranted. The total number of such star clusters in each galaxy is listed in Table 1. The positions of the identified star clusters in the plane-of-sky is shown in Fig. 1, overplotted on the HST image of the galaxy.

The magnitude limit of -6 in V that we set for visual inspection corresponds to a 1000 M , 6 Myr star cluster with colour excess E(B -V) = 0.25 (Calzetti et al. 2015a). Thus our catalogs are incomplete at lower masses. However, this limit corresponds approximately to the completeness limit of the underlying photometric catalog. Adamo et al. (2017) carry out artificial cluster tests for NGC 628 (distance of 10 Mpc), and find that LEGUS produces a complete cluster sample down to a cluster mass of 5000 M for cluster ages <200 Myr. Since our clustering results are driven by much younger clusters (<10 Myr), where we are complete down to even lower masses, incompleteness will have minimal impact on the results and analysis.

The final catalogs of which we make use in this work are publicly available online 2 on the Mikulski Archive for Space Telescopes (MAST) for all galaxies except NGC 3627 and NGC 5457. The catalog for these galaxies will be published in a forthcoming paper (Linden et al., in preparation).

It is important for our analysis of the spatial distribution of clusters to de-project the cluster positions from the plane of the sky to the plane of the galaxy. This is necessary, especially for higher inclination angles, as the inclination modifies the true spatial separations between the clusters in the plane of the galaxy. The first step in our analysis is therefore to de-project cluster positions.

To perform this correction, we assume that each galaxy can be described with an axisymmmetric flat rotating elliptical disc. We then correct the position of each star cluster in a two step process. First, we rotate the intrinsic positions of the clusters by an angle φ in the clockwise direction about the center of the galaxy, where φ is the position angle measured anticlockwise from the celestial north, to align the major axis of the galaxy in the north-south direction. For a cluster with RA and Dec. positions x and y relative to the center of the galaxy, with x increasing along RA (toward the left of Fig. 1) and y along DEC (toward north of Fig. 1), we compute the new positions as

where x and y are the position-angle corrected RA and Dec. of the cluster. Our second step is to correct for the line-of-sight inclination angle i by dividing x by cos i while leaving the y position unchanged.

Thus the final positions of all clusters are x i = x /cos i and y i = y . The values of φ and i that we use for each galaxy are provided in Table 1. The RA and Dec. of the galaxy centers are taken to be the reported sky coordinates for the galaxies in the NASA/IPAC Extragalactic Database (NED). 3 We use these de-projected cluster positions (x i , y i ) for all calculations of the correlation function in this paper.

To investigate the hierarchical distribution of young star clusters in galaxies, we use the angular TPCF 1 + ω(θ), where θ is the angular separation between a pair of clusters in the plane of sky. We refer to this quantity as the TPCF or the correlation function interchangeably for the remainder of the paper. The physical meaning of the TPCF is that, if one examines an annulus of radius θ and infinitesimal width dθ centered on one cluster, 1 + ω(θ) is the ratio of the probability of finding another cluster within this annulus to the probability of finding one in an identical annulus that is placed at a random position, rather than centered on a known cluster. Mathematically, we express this by writing the conditional probability dP(θ) that a pair of clusters in a region is separated by an angle θ as,

where N is the average surface density of clusters in the region per steradian, d 1 and d 2 are infinitesimal solid angle elements around clusters 1 and 2, respectively, and 1 + ω(θ ) has the typical form of a correlation function,

where N(θ 1 ) and N(θ 1 + θ ) are the local surface densities around two regions separated by angle θ , and the averaging is done over all such pairs of regions. From the definitions of the quantities above, we can interpret ω(θ) as the quantity that represents the excess probability above a purely random Poisson distribution of finding a pair of points separated by an angle θ . Indeed, for a purely random distribution, ω(θ) = 0, and by corollary, the correlation function 1 + ω(θ) = 1. For a clustered distribution ω(θ ) > 0, and for a scale-free clustered distribution such as a fractal, the correlation function 1 + ω(θ) is a pure power law with a negative slope, up to the scale at which the distribution remains a fractal (Calzetti, Giavalisco & Ruffini 1988).

A number of authors have proposed estimators to calculate ω(θ) for a given observed distribution of point-like objects (Peebles 1974;Peebles & Hauser 1974;Sharp 1979;Shanks et al. 1980;Hewett 1982;Hamilton 1992;Landy & Szalay 1993). We use the one proposed by Landy & Szalay (1993) as it attempts to correct for effects near the edge of the field of view, and provides estimates whose errors are largely Poisson distributed. It uses a combination of the data sample and a random sample that populates the field of view of the data. The Landy-Szalay estimator (LS, hereafter) is calculated as

where DD is the number of data-data pairs, DR is the number of cross-correlated data-random pairs, and RR the number of randomrandom pairs, counting all pairs in the range of separations θ ± dθ, where dθ is the adopted width of the discrete separation bin for which the TPCF is calculated. It is typically desirable to have a large enough random sample to cross-correlate with, so as not to introduce any additional Poisson error in the estimator. To accommodate this, we normalize the counted pairs at a separation θ to the total number of possible pairs in the data, random, and cross-correlated distributions. This leads to the following modifications to the definition of DD, DR, and RR

where N D and N R denote the total number of data and random points in the distribution. We ensure that the sky coverage and geometry of the random sample is as identical as possible to that of the data, to allow accurate TPCF computation, especially close to the edge of the field of view. This is done by preparing a Poisson distribution of points occupying the HST footprint of the observed galaxy, then masking out any unsampled regions (if any) that happen to fall on the chip gaps of the ACS instrument. The footprint of the galaxy is prepared from the observed V-band image using the FootprintFinderfoot_3 tool that is publicly available, and adjusted to account for the deprojection of the galaxy. We do not, however, mask potentially dust-extincted regions in the galaxy such as dust lanes, since Grasha et al. (2015) found that doing so does not affect the resulting TPCF.

To compute the value of the TPCF, we use our own modified version of the TPCF functionality offered by the PYTHON ASTROML 5 module (Vanderplas et al. 2012), which uses the scikit-learnfoot_5 library as a backend for fast computations of pairs at given separations (Pedregosa et al. 2011). We use a bootstrap method (Efron & Tibshirani 1994) with 100 bootstrap samples to estimate the value of and error bars on ω in each bin. The code we used to perform the analysis is publicly available on GITHUB. 7

At angular separations that approach the size of the telescope field of view, the TPCF is prone to bias due to edge effects. The Landy & Szalay (1993) estimator attempts to compensate for this by smoothing the steep fall in the TPCF expected near the edge, as the number of pairs in the data goes to zero near the boundary. However, this correction is not perfect, and hence separations where it is significant should be interpreted with caution. To obtain an estimate for the scale where this correction starts to matter, we perform numerical simulations of toy 2D fractal distributions truncated by square fields of view of different sizes, and estimate the size l edge at which edge effects become significant. We provide details of the procedure we use to estimate this scale in Appendix A. We find that l edge ∼ R max /5, where R max is the size of the field of view, and the value of the TPCF for scales beyond l edge has a significant contribution from smoothing by the estimator we use. Thus, in all our plots of TPCFs, we use grey shading to indicate separations θ > θ max /5, where θ max is the angular extent of the deprojected HST field-of-view footprint. In many cases. the de-projected field of view is not a square, in which case we use the longest side for θ max . The separation θ max /5 is used to caution the reader about scales where edge effects might play a role and the TPCF should be interpreted with caution.

We compute the angular TPCF 1 + ω(θ) for star clusters in the galaxies of our sample through the method outlined in Section 3, using 20 angular separation bins that are logarithmically spaced in the range 10-5000 pc. We use the same bins (in linear rather than angular units) across all the galaxies to allow for consistent comparison of the TPCFs. Our choice to use 20 bins represents a compromise between resolving finer features of the TPCF and avoiding excessive shot noise. We mask bins for which the value of 1 + ω ≤ 0, which we sometimes encounter at large separations when the number of pairs in the data are very low. We also do not show the TPCFs for bins where the median value of ω is lower than the bootstrap-calculated error in ω, a condition that arises occasionally for narrow, low separation bins where the number of pairs is small.

In order to isolate evolutionary effects, for the bulk of this paper, we will consider only TPCFs for clusters separated by age. Separating clusters by age is important because there are also physical differences between the galaxies that lead to perceptible variations in their TPCFs, even within the same age group. Computing the TPCFs of clusters of all ages as a single distribution tends to mix the two causes of variation, making it impossible to disentangle variations in galaxies spatial structure from variations in their star formation histories. However, for reader convenience, we do present the TPCF for the combined cluster sample without binning by age in Appendix B. For the remainder of our analysis, we choose an age of 10 Myr to separate young (T 10 Myr) and old (T > 10 Myr) clusters. This choice is motivated by recent evidence that starforming giant molecular clouds are dispersed by feedback from massive stars in 1-5 Myr (see, Chevance et al. 2020, and references therein), where we have chosen a conservatively higher value as some galaxies have a low number of clusters at ages T < 5 Myr, and our age estimates carry some uncertainty. Thus, our division of young and old clusters corresponds roughly to those that are probably still associated with their natal molecular cloud, and those that are not, respectively. While in principle this time-scale would be different in each galaxy, we choose a consistent value for simplicity. In Fig. 2, we show the resulting TPCFs for young and old clusters in the galaxies. We note that NGC 3344, NGC 5253, and NGC 7793 have very few clusters with ages T > 10 Myr, and thus their TPCFs are extremely noisy even for the bins that have non-zero correlation. For this reason, we omit the TPCF for older clusters in these galaxies from the figure.

We find that, in general, there are significant variations between the TPCF of young and old clusters for a given galaxy. For instance, younger clusters seem to have relatively higher values of correlation and seem to show, at least qualitatively, power-law behaviour (straight line in a log-log plot) as a function of separation θ over a range of scales. Older clusters, on the other hand, seem to show TPCFs that are relatively flat at small separations and show some form of smooth fall-off at large separations, suggesting some sort of evolutionary effect. However, there is also significant variation in the TPCFs for a given age group among the galaxies, especially for the young clusters group. While some galaxies show powerlaw behaviour in their young cluster TPCFs over the entire range of scales where we measure it, others seem to show a steep power law at small scales, followed by a shallow/flat power law at larger scales, and a relatively sharp break between the two regimes. In addition, there are some galaxies, specifically the dwarfs, that seem to show power-law behaviour at small separations, followed by a smooth fall-off at larger separations. These differences, we expect, should be due to physical differences in the underlying star cluster distribution, which is presumably set by the host galaxy. We intend to understand both forms of differences, evolutionary and physical, in the sections below. First, in Section 4.2, we formally characterize the qualitative features we described above in the observed TPCFs, by fitting functional forms that reproduce these features, and using model comparison to choose the form that best describes the TPCF. We then attempt to infer the underlying distributions that might give rise to the observed TPCFs through the use of physically motivated toy model distributions. This allows us to infer the changes in the underlying star cluster distribution with host galaxy and with age. Following that, in Section 4.3, we probe evolutionary effects in further detail for a single galaxy, by adding more age groups, and measuring the TPCFs for clusters that fall in them, in order to obtain finer time resolution.

In this section, we attempt to explain the observed features in the TPCFs outlined in Section 4.1 by characterizing them quantitatively (Section 4.2.1), and comparing with physically motivated toy model distributions to infer the underlying star cluster spatial distributions. The full description of the parameters of the toy models and how they influence the TPCF of the distribution is provided in Appendix C. Here, we just briefly describe and motivate the toy models for each identified feature, and infer the properties of the spatial distribution of the clusters from the comparison between the toy models and our measured TPCFs.

The first step to understanding the TPCFs is to classify them by morphology of their features. We do this by defining three functional forms/models for the TPCF that represent qualitatively the three features described above, attempting to fit them to the observed TPCFs, and performing a statistical comparison of the fits to identify which functional model describes the data best. The three models can be described qualitatively as, (i) Model S: A single power law with a fixed slope. (ii) Model PW: A piecewise power law consisting of two fixedslope segments separated by a transition/break point.

(iii) Model PF: A power law with an exponential cutoff to represent a power law that falls off smoothly at larger scales.

Quantitatively, we define model S as

where A 1 is the amplitude and α 1 the power-law slope. Model PW is given by

where α 1 and α 2 are the two power-law slopes, β the transition point, and A 1 and A 2 the amplitudes, which are related by A 2 = A 1 β α 1 -α 2 to ensure continuity at θ = β. Lastly, Model PF is given by

where the slope is as defined for the single power-law case, and θ c , is the scale above which the TPCF falls off.

We fit all three model TPCFs to the observed correlation function 1 + ω(θ) of every galaxy in our sample using a Markov Chain Monte Carlo (MCMC) method (MacKay, Kay & Press 2003). The parameter vectors to fit for are λ S = (A 1 , α 1 ), λ PW = (A 1 , α 1 , α 2 ), and λ PF = (A 1 , α 1 , θ c ) for Model S, PW, and PF, respectively, and the likelihood function is given by

where D(θ) is the observed value of 1 + ω(θ) at separation θ , M(θ |λ) the corresponding model value at this separation for the parameter vector θ and σ D (θ) the error in the observed TPCF value. The priors we use on the parameters are: A 1 > 0, -5 < α 1 ≤ 0, -5 < α 2 ≤ 0, θ min < β < θ max , and θ min < θ c ≤ 5θ max , where θ min and θ max are the minimum and maximum separations over which we compute the TPCF. We use the PYTHON package EMCEE (Foreman-Mackey et al. 2013) to perform the calculation, using 300 walkers, with a total of 5000 steps, discarding the first 200 steps as burn-in. We verified that the MCMC had reasonably converged in such a case through visual inspection of the MCMC chain.

To determine which model among the three described above is the best description for a given galaxy and age group, we compute 10 Myr (blue) and T > 10 Myr (red) for each galaxy, with the best-fitting functional form model for both age groups overplotted (dashed lines) in their corresponding colours. The functional models that are fitted are: Single power law (Model S), PieceWise power law (Model PW), and power law with an exponential fall-off (Model PF), and are represented as dotted, dashed, and dot-dash line styles, respectively. The functional form among these that fits the TPCF best for both age groups is reported in the legend, with the best-fitting parameters and superior model reported in Table 2, and the detailed fitting procedure outlined in Section 4.2.1. The bottom x-axis denotes θ , the angular separation in arcsec, and the top axis denotes δx, the corresponding linear separation in parsec (pc) using the distance to the galaxy reported in Table 1. Grey shaded regions denote the estimated range of separations where edge effects could play a role in the TPCF (see Section 3.3). 7), PieceWise power law (PW) -equation ( 8), and Single power law with exponential Fall-off (PF) -equation ( 9). The parameters include: α 1 , the power-law slope common to all models, α 2 the slope of the second power law in Model PW, β the transition point in arcsec ( ) for Model PW, and θ c the scale separation in arcsec of the exponential fall-off in Model PF. The AIC values and the best-fitting parameters are obtained from MCMC fits of the three aforementioned models to the observed TPCFs. Further description of the models, their parameters, and the fitting procedure are given in Section 4.2.1. The AIC is given by equation ( 11) and the model with the lowest AIC value is considered the best-fitting model. This best-fitting model is reported in column 6 ('Best Model') and only the best-fitting parameters associated with that model are reported. Galaxies that had very few old clusters (i.e. NGC 3344, NGC 5253, and NGC 7793), and as a result very noisy TPCFs do not have AIC values or fits reported. In addition, in some cases, a model has a very low likelihood, which leads to extremely high AIC values, in which case we denote the AIC value with a lower limit of 10 3 , i.e. >10 3 .

the Akaike Information Criterion (AIC; Akaike 1974) for each fitted model. The AIC is an estimator for the relative quality of statistical models, given a set of data, which compares goodness of fit along with a factor that penalizes for a higher number of parameters, preventing over-fitting a model to data. Given a set of candidate models, with their respective AIC values, the preferred model is the one with the minimum AIC value. In this study, we use the so-called corrected AIC (Hurvich & Tsai 1989), which adds a correction term to the traditional AIC value, making it suitable for small sample sizes, and converges to the traditional AIC for an infinite sample size (Burnham & Anderson 2004). It is given by

where N λ is the number of parameters for a model, and N θ is the number of angular separation bins for which the correlation function is calculated. In Table 2, we report for each galaxy and age group (i.e.

young and old) the AIC values for all three models, the model with the minimum (best) AIC, and the resulting best-fitting parameters for the superior model. In Fig. 2, we overplot this best-fitting model for both young and old clusters on their corresponding measured TPCFs with different line styles, and list the best-fitting model name in the legend for each galaxy and age group. We find that the best-fitting model for the older clusters in a galaxy is typically Model PF, and this is more or less consistent across the sample. On the other hand, the best-fitting model for younger clusters varies. Among the spirals, the TPCFs of NGC 1313, NGC 1566, NGC 3627, and NGC 5194 are best fit by Model S, whereas NGC 0628, NGC 3344, NGC 5457, NGC 6503, and NGC 7793 prefer Model PW. For the dwarf galaxies, NGC 3738 and NGC 5253 prefer Model PF for both young and old clusters, and NGC 4449 is best fit by Model PW and Model PF for young and old clusters, respectively. To understand these differences and their physical implications, it is important to first identify what sort of underlying star cluster distribution gives rise to the three models. We perform this exercise for all three fit models by comparing with physically motivated toy model distributions below.

Many previous studies have reported pure power law TPCFs of the form described by model S, not just for star clusters, but also for individual stars and for tracers of gas (Zhang et al. 2001;Bastian et al. 2005;Scheepmaker et al. 2009;Gouliermis et al. 2014Gouliermis et al. , 2015bGouliermis et al. , 2017;;Grasha et al. 2015Grasha et al. , 2017aGrasha et al. , 2018Grasha et al. , 2019;;Shabani et al. 2018). We interpret power laws in the TPCFs as a sign of a self-similar hierarchical distribution in the star clusters. This is because the TPCF of a self-similar fractal distribution of points is a pure power law of the form 1 + ω(θ) ∝ θ α (Calzetti et al. 1989;Larson 1995), with the power-law slope α related to the 2D fractal dimension D 2 as D 2 = α + 2. In Appendix C1, we verify that toy fractal distributions show a pure power law TPCF up to the largest scale of the hierarchical structure, beyond which the TPCF flattens to approach a value of 1 + ω(θ) ∼ 1. The amplitude A 1 of the power law (see equation 7) depends on both D 2 and the field of view length scale R s , and agrees reasonably well with the analytical predictions of Calzetti et al. (1988). Physically, a self-similar hierarchy in the star clusters can originate from fractal density distributions in the natal star-forming gas, which in turn result from supersonic turbulent motions in the ISM (Elmegreen 1993;Mac Low & Klessen 2004;Kritsuk et al. 2007;Federrath et al. 2009).

The young clusters in some spiral galaxies show a piecewise power law TPCF. By examining the fits for these galaxies in Table 2, we can clearly see that the best-fitting slope is steeper at smaller scales, and shallower at larger scales, i.e. α 1 < α 2 (except in NGC 4449). The angular separation β where the slope changes lies in the range ∼2-14 arcsec, which corresponds to scales of ∼80-500 pc. The break in slope cannot be edge effects, since it occurs at a scale well below the range at which we expect edge effects to play any role. Another possible cause for a break in the TPCF is a 2D-to-3D transition of the underlying distribution at the scale height of the galaxy. In other words, the distribution of star clusters could be 3D at separations smaller than the scale height, but at separations beyond the scale height any pair of clusters would both lie on the plane of the galaxy, rendering the distribution 2D. The projected fractal dimension, which is related to the slope of the TPCF, is different if we are looking at the projection of a thin slice (2D) rather than a thick disc (3D), and this projection effect leads to a change in the slope (see, Sánchez & Alfaro 2008, for detailed models). Observations of neutral hydrogen, far-infrared dust, and γ -ray emission in nearby external galaxies are consistent with this mechanism [Elmegreen, Kim & Staveley-Smith 2001;Miville-Deschênes et al. 2003;Ingalls et al. 2004;Dutta et al. 2009;Szotkowski et al. 2019;Besserglik & Goldman 2021, although see Koch et al. (2020) for an alternative explanation for this transition]. However, our cluster TPCFs are not: as discussed in Sánchez et al. (2010), this transition should lead to a steeper slope at larger separations and a shallower ones at smaller separations, which is the opposite of what we find.

Having ruled out edge effects and scale height effects, we conjecture that the breaks we see in galaxies with PW-type TPCFs represent real transitions from a fractal distribution at smaller scales set by turbulence to a mostly random distribution at larger scales where 2D galactic dynamics become more important than turbulence. Such a transition produces a shallow slope at large separation and a steeper slope at small separation, which is what we observe, and what is also seen in the correlation function of stars and H II regions in M33 (Odekon 2008;Sánchez et al. 2010). In Section C1, we test this scenario by creating toy fractal distributions that are scale-free up to a maximum scale L max , and which are Poissonian at larger scales. We vary L max in our models, attempt Model PW fits to them, and find that the transition point parameter β picks out the randomization scale L max quite well. Hence we infer Model PW fits to represent fractal distributions that are scale-free only up to some maximum size scale ∼β, and become non-fractal at larger-scales. The match between the toy models and observations is not perfect, however: the measured TPCFs do not transition sharply to a completely flat slope like the toy models, but rather more smoothly to a value of 1 + ω ∼ 1. We speculate that this is because the distribution beyond the transition point β is not entirely Poissonian, since the large-scale distribution of clusters in a galaxy is clearly non-uniform on scales approaching the galactic scale length. Indeed, in the following section, we find direct evidence for this effect.

The third class of model, i.e. a power law that smoothly transitions to an exponential, is the best fit for the old clusters in all galaxies where we have enough old clusters to carry out a fit, and is also the best fit for young clusters in some of the dwarf galaxies in our sample. Following Mao et al. (2015), we hypothesize that this functional form reflects the large-scale distribution of the clusters in a galaxy, to which the young clusters converge as they age. To test this hypothesis, we distribute clusters in our toy models using a standard large-scale distribution: a radially thin exponential disc with a given scale length, and a Gaussian distribution in the vertical z direction with a characteristic scale height. While there are more detailed models to describe the large-scale distribution of clusters in a galaxy, we chose the exponential disc model for simplicity. The toy model and its parameters are described in Appendix C2. We find that the TPCFs of the toy models display a smooth fall-off with separation, with a Model PF fit yielding θ c that corresponds to the exponential scale radius r c of the radial distribution (see, right-hand panel of Fig. C3). Adding logarithmic spiral arms in the azimuthal direction to the exponential disc does not significantly change the TPCF. Since there is a close resemblance between the observed TPCFs and these toy models, we conjecture that PF-type TPCFs simply reflect the large-scale radial distribution in galaxies, which is reasonably well-described by a thin exponential disc with radial scale length r c approximately equal to the fitted scale length θ c .

We have shown that the TPCFs for young (<10 Myr) and old (>10 Myr) clusters are significantly different. In most galaxies, older clusters have flatter, lower amplitude TPCFs, indicating that both they are less clustered overall, and that, unlike young clusters, older clusters do not follow a scale-free fractal structure. While a number of authors have reported qualitatively similar results, (Odekon 2006;Sánchez & Alfaro 2009, 2010;Grasha et al. 2015Grasha et al. , 2017aGrasha et al. , 2018Grasha et al. , 2019)), limited sample sizes have made it difficult to follow the evolution of the TPCF with time in detail. Because our sample of star clusters is among the largest available for this type of analysis, we can, at least for some galaxies, bin by age much more finely, and thereby obtain a higher resolution picture The bestfitting model for each age group is overplotted with dashed lines, and the model name and parameters for each curve reported in the legend, both coloured by age group. It is evident that the TPCFs of younger clusters are scale-free power laws that decrease in slope as the clusters age. Older clusters are distributed more evenly across the disc and hence show shallower powerlaw slopes at small scales, with an exponential fall-off at large scales due to the correlation imposed by the overall radial scale in the disc. of TPCF evolution. We therefore divide star clusters into four age brackets: T < 2 Myr, 2 < T < 10 Myr, 10 < T < 100 Myr, and T > 100 Myr, and compute the TPCF for the distribution of clusters in these age brackets. We then use the fitting method outlined in Section 4.2.1 to choose the best-fitting functional form that describes the TPCF. We can only carry out this analysis for a subset of galaxies in our sample, as the rest do not have enough clusters at a wide enough range of ages.

We show the result of this analysis in Fig. 3 for NGC 5194, which has the highest number of clusters among all the galaxies in this study. We find similar qualitative behaviour for NGC 1313, NGC 1566, and NGC 0628, the other galaxies in our sample for which we were able to perform this analysis, albeit with substantially larger uncertainties due to the smaller numbers of clusters available. In the case of NGC 5194, we find that the best-fitting model for clusters in the younger two categories are clearly single power laws (Model S) with a slope α 1 that decreases from -0.55 for the youngest clusters (T < 2 Myr) to -0.38 for clusters with 2 < T < 10 Myr. We do not see any signs of an exponential fall-off from the exponential disc distribution for the younger clusters, or a break in the power law that might indicate an outer limit to the scale-free structure. This is consistent with the visual impression from Fig. 1, which shows that younger clusters are mostly concentrated in hierarchically structured patterns that predominantly seem to trace the spiral arms in spiral galaxies, or the central regions of dwarf galaxies. On the other hand, clusters with ages 10 < T < 100 Myr are best fit by a power law with an exponential fall-off at large scales (Model PF), and a shallow power-law slope of α 1 = -0.28 at small scales. This marks a transition phase where clusters are losing their natal fractal structure and thus have a shallower power law. However, these clusters are also old enough that they are distributed fairly uniformly across the extent of the disc, such that the imprint of the overall exponential radial distribution becomes evident at larger Figure Schematic summarizing the three functional forms fitted to the TPCF and the physical quantities that can be inferred from them, obtained based on our analysis using toy model distributions outlined in Section 4.2. The three functional forms fitted are: Single power law (Model S: equation 7), Piecewise power law (Model PW: equation 8), and Power law with exponential fall-off (Model PF: equation 9). α 1 is the small length-scale power-law slope in all three models, β is the transition scale in Model PW, θ c is the exponential scale separation of Model PF, and θ max is the largest bin to which the TPCF is measured. From these models and their best-fit parameters, we infer values of l corr , the largest scale to which hierarchical structure extends (see Section 4.4.1), D 2 , the 2D fractal dimension of the fractal distribution (see Section 4.4.2) -both calculated from the young cluster TPCF -and r c , the exponential scale radii of the large-scale distribution in the galaxy, which is calculated from fits to the old clusters TPCF. Our estimate for l corr is taken to be ≈β from Model PW, and θ max from Model S. Our estimate of D 2 = α 1 + 2 for all three models, and that for r c is estimated to be ≈θ c . scales. Finally, the oldest clusters in the galaxy (T > 100 Myr) have a negligible small-scale slope α 1 ≈ 0, and a clear exponential fall-off at large scales, suggesting that for this age group fractal structure at all scales is completely lost, and the only remaining contribution to the TPCF comes from the large-scale radial structure of the disc.

This result is consistent with the finding in Grasha et al. ( 2019) that star clusters in NGC 5194 become spatially decorrelated from molecular clouds by ages of ∼50-100 Myr. While our results on the TPCF in other galaxies are too noisy for us to perform a similar measurement in them, we note that Grasha et al. (2018) found a lower cluster-molecular cloud decorrelation time in NGC 7793. Thus it is likely that the cluster-cluster decorrelation time that we are measuring will also depend on the host galaxy and its environment.

The three functional forms that we find provide a good description of the cluster TPCFs -models S, PW, and PF -and are characterized by three parameters: the largest scale up to which there is fractal signatures in the distribution l corr , the 2D fractal dimension of the distribution D 2 in the range of separations up to l corr , and the scale radius r c beyond which the TPCF declines exponentially. We provide a schematic summary of these quantities, and their relationship to our functional forms, in Fig. 4. In the remainder of this section, we investigate the distribution and variation of each of these quantities over the galaxy sample, and discuss possible physical origins for their values. -See Fig. 4 for a schematic outlining the method we use to obtain the values above from the fits in Table 2. D 2 : 2D fractal dimension inferred from the TPCF of young clusters, with error bars obtained from the fits. l corr : Largest scale of hierarchical structure inferred from the young cluster TPCF. Error bars for l corr , if any, take into account the uncertainty in the distance to the galaxy and the fit uncertainty. The estimates of l corr for galaxies where the TPCF is best fit by Model S are lower limits. We do not report a value of l corr for NGC 3738 since it shows no evidence of fractal structure (D 2 ∼ 2). r c : Exponential scale radii derived from the TPCF of older clusters. Error bars take into account distance and fit uncertainties. We do not report r c values for NGC 3344 and NGC 7793, due to their lack of older clusters. For the same reason, for NGC 5253, r c is obtained from the Model PF fit to the young cluster TPCF.

The largest scale of hierarchical structure, l corr , denotes the maximum separation up to which star clusters are distributed in a scale-free fractal distribution, and beyond which star clusters are uncorrelated with each other. Since star clusters form from the underlying gas distribution in the ISM, we expect that l corr is also approximately the size of the largest coherent gas structures (Efremov 1995). We infer l corr from the TPCFs of the young clusters (T 10 Myr), since these have been least influenced by evolutionary effects, and thus should most closely reflect the distribution at cluster formation. The method we use to estimate this scale for a galaxy depends on the best-fitting model for its young cluster TPCF (see Table 2), and is summarized in the schematic shown in Fig. 4. For Model PW galaxies, we take l corr to be the TPCF transition point β beyond which the distribution of clusters randomizes. For Model S galaxies, where the power law extends out to the last bin of measurement, we can only estimate a lower limit for this scale, taken as the largest scale for which we can measure the TPCF. For the dwarf galaxies, namely NGC 4449 and NGC 5253, we estimate l corr as the scale where the power law sharply turns down due to the effect of the exponential disc distribution at larger scales (see Fig. 2); this estimate is likely an upper estimate, as the effect of the exponential disc could be present even at smaller scales, and it difficult to disentangle the power law part from the exponentially falling part of the TPCF. We do not calculate a value of l corr for NGC 3738 since it does not show any sign of scale-free fractal structure at any scale, since α 1 is found to be ∼0. The values of l corr are reported in Table 3.

As we can see, the values of l corr vary among the galaxies and lie in a rather broad range from ∼ 100 pc in NGC 7793 to upwards of 3000 pc in NGC 5194. We compare these values with the various galaxy properties listed in Table 1: the standard isophotal radius R 25 , the morphological T type, the galaxy stellar mass M * , the UV-derived star formation rate SFR UV , and the stellar mass and star formation rate per unit area ( * and SFR ). We show scatter plots of l corr against these quantities in Fig. 5, and report the Pearson correlation coefficients ρ and their corresponding p -values in the Figure legend. The value of ρ for a pair of variables lies in the range -1 to 1, with 1 (-1) indicating perfect linear correlation (anticorrelation) and 0 denoting no linear correlation; p is the probability of obtaining a correlation coefficient ≥ρ from a pair of variables that have, in fact, zero correlation (i.e. the null hypothesis), and is thus a measure of the statistical significance of the measured correlation. A value of p < 0.05, meaning < 5 per cent probability of a false positive, is typically interpreted as statistically significant (see, for instance, Freedman, Pisani & Purves 2007). We find moderately significant correlations of l corr with M * (ρ = 0.65, p = 0.03), SFR UV (ρ = 0.56, p = 0.07), and SFR (ρ = 0.69, p = 0.02), and no significant correlation with other quantities. The three detected correlations strengthen if one considers only the spirals in the sample. This analysis suggests that more massive and brighter galaxies (which also have higher star formation rates) tend to contain correlated complexes undergoing hierarchical star formation with larger sizes than are found in less massive galaxies, and agrees with similar signs of correlation found using the TPCF of star clusters in Grasha et al. (2017a). It is interesting to note that galaxy size (R 25 ) shows no correlation with l corr (ρ = 0.17, p = 0.62), and that the areaaveraged star formation rate ( SFR ) correlates more strongly with l corr than the total star formation rate SFR UV . This suggests that l corr is determined more by the physical conditions of the star-forming gas than by the overall size of the galaxy. However, we note that ρ and the associated p-values are calculated using the lower (upper) limit value in the case of Model S (PF) fits, and hence might be different if we had real values. In addition, we caution that our sample is limited to only 12 galaxies, so any correlations are only suggestive, not conclusive, due to the low sample size.

What physical mechanisms set l corr ? The correlation with M * and SFR suggests that the gravitational potential of the matter (stars and gas) in the galaxy is important, with stronger potentials leading to larger complexes. This agrees with the physical picture of starforming clouds at galactic scales being formed through gravitational instability in the disc, with their size ultimately limited by some topdown mechanism that prevents them from growing too large. Galactic rotation is a candidate stabilizing mechanism on large scales, and such a picture would be consistent with the results of Grasha et al. (2017b), who find that a velocity gradient set by shear could explain the variation among the galaxies in the largest scale up to which pairs of star clusters are correlated in age. In this scenario, gravitational instability is unable to create structures past a certain maximum size, beyond which galactic rotation stabilizes the disc. The natural scale in this case is the Toomre length l T (Toomre 1964;Escala & Larson 2008),

where G is the gravitational constant, g is the gas surface density, and κ the epicyclic frequency of rotation. For a flat rotation curve, which we assume here, κ = √ 2 , and = v rot /r is the angular rotational velocity calculated from the flat rotational velocity v rot at a given galactocentric radius r. To check this hypothesis, we compute l T for all the spiral galaxies in our sample; we omit dwarf/irregular galaxies due to the lack of robust observed rotational curves as it 3 with the isophotal radius R 25 , stellar mass M * , UV-derived star formation rate SFR UV , morphological T-value, stellar mass surface density * , and star formation rate surface density SFR of the host galaxy. The different marker styles denote the three ways that l corr is estimated (see Section 4.4.1 and the schematic in Fig. 4) based on which functional form fits the young cluster TPCF best (reported in Table 2). The Pearson correlation coefficient ρ and corresponding p-values of the correlation are provided for each pair of variables. We find signs of correlation of l corr with M * , SFR UV , and SFR , with stronger correlation if we restrict the sample to the spirals only. However, note the caveat that we use lower limits in the case of Model S galaxies to calculate the values of ρ and the associated p-values, and they may be different if we had constrained values of l corr instead of lower limits.

is unclear to what extent the dwarfs have a disc-like structure. This calculation requires estimates for g and v rot , and an appropriate choice for r. To calculate g , we use galaxy-averaged molecular gas (H 2 ) surface densities reported in the literature where available, and estimated total molecular gas masses M H 2 from the literature divided by πR 2 25 otherwise. We use H I rotation curves and their reported rotational velocities available in the literature to infer v rot . We choose the representative radius r at which to calculate l T to be the median galactocentric radius of the young star clusters in our star cluster catalogs. We show l T versus l corr in Fig. 6. We find a reasonably strong (ρ = 0.75) and statistically significant (p = 0.01) correlation between the two. The two scales we calculate, although correlated, are not identical; in general l corr is larger than l T by a factor of a few. We caution that modern treatments of the Toomre instability include the effects of multiple stellar populations along with the gas, the effects of finite thickness, and the dissipative nature of gas (see e.g. Romeo & Falstad 2013). However, we lack measurements of the stellar velocity dispersion or disc scale height, which would be required to include these effects, and thus we limit our comparison to the simple pure-gas Toomre length. It is important to extend this comparison to a larger sample of galaxies, to obtain more robust and conclusive results for the importance of such a mechanism.

Another quantity of interest is the fractal dimension of the hierarchical distribution at scales l < l corr . The fractal dimension is a quantity 12), and the inferred largest scale of hierarchical structure in star clusters l corr for the spiral galaxies in our sample. Marker styles are as outlined in Fig. 5. We find a statistically significant correlation with a Pearson correlation coefficient ρ = 0.75 and a p-value of 0.01. This qualitatively suggests a physical picture where the largest scale of the hierarchy in star clusters is set by rotation-supported gravitational instability of the gas lying in the galactic disc. 5, but for the inferred 2D fractal dimension D 2 of the young cluster distribution, except for NGC 3738 that has a value of D 2 corresponding to a Poissonian distribution. Note that some points have very small errors, which are not visible. Overall, we find weak signs of correlations for D 2 , which are, however, not statistically significant. that characterizes self-similar structure in a distribution, with lower values corresponding to less space-filling hierarchical structures. Self-similar hierarchies are proposed to be set self-consistently by interstellar turbulence in the ISM gas (Elmegreen & Scalo 2004;Federrath et al. 2009). If this picture is correct, the result should be a nearly universal value for the fractal dimension, as has been proposed in earlier studies (Feitzinger & Galinski 1987;Elmegreen & Falgarone 1996). Previous studies of the fractal dimension of the gas and/or dust distribution in galaxies have generally been consistent with the hypothesis of a universal fractal dimension (see, Table 1, Shadmehri & Elmegreen 2011). However, Sánchez & Alfaro (2008) find statistically significant variation in the fractal dimension of H II regions with host galaxy and/or environment. Here, we investigate whether the fractal dimensions of the distributions of star clusters in our sample are the same in all galaxies, and if not, how its variation correlates with other galactic properties.

We compute the 2D fractal dimension D 2 from the power-law slope of the fits to the young cluster TPCFs reported in Table 2, using the relation D 2 = 2 + α 1 , where α 1 is the fitted slope of the power law. Since all three fit models include α 1 as a parameter, we obtain a corresponding D 2 for all galaxies in our sample. This approach is summarized in the schematic shown in Fig. 4. As in the previous section, we do this for the young clusters TPCF, which should more closely reflect the fractal dimension of the natal gas supposedly set by interstellar turbulence. We list D 2 for each galaxy in our sample in Table 3. We find variations well beyond the computed 1σ errors, with D 2 lying in the general range 0.5-1.6, with the exception of NGC 3738, which has a value D 2 corresponding to a completely random distribution, i.e. D 2 ∼ 2.0. This suggests that, consistent with Sánchez & Alfaro (2008), and contrary to earlier suggestions (Feitzinger & Galinski 1987), the hierarchical structuring in the star cluster distribution does not show signs of universality and depends on the host galaxy and its properties in a way that the gas distribution apparently does not (Elmegreen & Falgarone 1996;Shadmehri & Elmegreen 2011).

We show scatter plots of D 2 versus various galaxy properties in Fig. 7; we report the Pearson correlation coefficient for each of the comparisons shown in the corresponding figure panels. As with l corr , we find at most marginal evidence for correlation of D 2 with M * , SFR UV , and SFR ; the correlation is stronger if we consider only the spirals in the sample, but remains below the level of statistical significance. To the extent that we interpret the vague hints in our data, they suggest that more massive galaxies have larger fractal dimensions (more space-filling distributions) than less massive galaxies. Such a trend for the inferred fractal dimension have been reported in earlier studies -i.e. brighter galaxies -quantified by their Bband absolute magnitude -have higher fractal dimensions than fainter ones (Parodi & Binggeli 2003;Odekon 2006;Sánchez & Alfaro 2008). In addition, Sánchez & Alfaro (2008) found that this correlation disappears when the irregular galaxies are included in their analyses, as irregular galaxies have fractal dimensions similar to the brightest spiral galaxies, but are also significantly fainter then them, qualitatively similar to what we find. It would be interesting to search for a similar effect for clusters using a larger sample of galaxies.

We also point out that there are earlier estimates for D 2 in the literature for a few of our galaxies. The values we obtain are consistent within the uncertainty in some galaxies, but not for all. For instance, Scheepmaker et al. (2009) estimate D 2 ∼ 1.6 for clusters younger than ∼ 30 Myr in NGC 5194, which is consistent with our result (1.6 ± 0.1). On the other hand, the values we obtain for NGC 0628 (∼0.9) and NGC 6503 (∼1.4) are different than earlier values 3) and r Spitzer , the value reported in the Spitzer Survey of Stellar Structure in Galaxies (S 4 G, Salo et al. 2015). Error bars are plotted for our inferred value r c , taking into account errors from the fit to the TPCF and the uncertainty in the distance to the galaxy. The error bars for r Spitzer only take into account the uncertainty in the distance as Salo et al. (2015) do not report error bars for their calculated scale lengths. A one-to-one dashed line (purple) is added to guide the eye. As we can see, r c reasonably reproduces r Spitzer for smaller galaxies, but overestimates it for larger galaxies where the HST field of view does not adequately cover the outer galaxy (see main text).

quoted for them in literature -i.e. 1.5 for NGC 0628 (Elmegreen et al. 2006;Gusev 2014) and 1.7 for NGC 6503 (Gouliermis et al. 2015b). This difference could occur for several reasons. For instance, these studies do not look at the hierarchical structuring of star clusters, but rather star-forming regions (in NGC 0628) or young stars (in NGC 6503), and there is no reason to assume that these structures all have the same fractal dimension. In addition, the NGC 0628 studies inferred a value of D 2 from the slope of the cumulative size distribution of star-forming regions, whereas we infer D 2 from the slope of the TPCF of young star clusters, a very different method. It is also well known that differential clustering estimates -such as the TPCF -are well suited to determining scales at which a change in clustering strength takes place (see for instance, Sánchez & Alfaro 2010). This feature, combined with our Bayesian approach to fitting various functional forms and hence slopes is important, especially in the cases of NGC 0628 and NGC 6503, which were best fitting by Model PW, and for which a fit to Model S only (analogous to the procedures used in earlier work, which implicitly assume a single power-law correlation function) would yield a significantly shallower slope, and hence a higher D 2 .

The exponential scale angle θ c , corresponding to a linear distance r c , is set by the radial distribution of clusters in the galaxy, and obtained by fitting Model PW to the TPCF of old clusters (T > 10 Myr), as indicated in the schematic shown in Fig. 4. We report values of r c in Table 3; note that the reported uncertainties include the uncertainty in the distance to the galaxy. We compare our r c values with r Spitzer -the scale radius of the galaxies in our sample estimated with the 3.6 and 4.5μm Spitzer Survey of Stellar Structure in Galaxies (S 4 G, Salo et al. 2015) in Fig. 8. We find reasonable agreement for galaxies that have lower values of r c , especially the dwarfs, but for most larger galaxies, we find r c r Spitzer . Why might this be the case? One possibility is that there are substantial uncertainties in r Spitzer , since S 4 G provides no estimate of uncertainties apart from those arising from the distance uncertainty. However, this seems unlikely to account for the factor of 3-4 discrepancy we find for large r c . A more likely explanation is that the HST field of view does not encompass the entire extent of the disc as it does for the smaller galaxies. To test whether this could lead to overestimates of r c , we artificially place a limited field of view on our toy galaxy models (see Appendix C2). We then compute and fit model PF to the TPCFs, and check whether the value of r c derived from the fitted θ c overestimates the true input value of the scale length we provide. In Fig. C4, we show that this is indeed the case: limiting the field of view to two galactic scale lengths leads to an overestimate of r c by a factor ∼3, roughly the observed discrepancy. We therefore tentatively conclude that the exponential cutoff found in Model PF gives a reasonable estimate of the scale length of the host galaxy, but only as long as the footprint within which the clusters are sampled extends to sizes significantly larger than the galactic scale length.

In this study, we investigate the hierarchical spatial distribution of young star clusters in 12 local galaxies cataloged with the LEGUS survey (Calzetti al. 2015a), using the angular TPCF 1 + ω(θ) as a function of angular separation θ. Our sample consists of various types, from irregular dwarfs to grand design spirals, allowing us to probe the effects of the host galaxy environment on the star cluster distribution. Estimated ages for the clusters obtained as part of the survey also allow us to study how the cluster distribution changes with age. We show that the TPCFs in all our galaxies are reasonably well-described by a model characterized by three parameters: the largest scale of hierarchical structure l corr , the 2D fractal dimension of the young star cluster distribution D 2 , and the radial exponential scale radii r c of the star clusters. We study how these parameters vary with the properties of the galaxies, to investigate the physical mechanisms that might be responsible in setting them. Our main results are summarized below.

(i) The TPCFs of younger clusters show large correlation amplitudes and strong fractal structure characterized by scale-free powerlaw TPCFs for separations θ ࣠ l corr . The TPCFs of older clusters, on the other hand, show shallow power laws, characteristic of more randomized distributions at smaller separations, and an exponential fall-off at larger separations (Fig. 2). Comparison with toy models shows that this fall-off is consistent with the cluster distribution following an overall exponential decline with galactocentric radius (Section 4.2.4).

(ii) The star cluster distribution loses its natal hierarchical structure gradually with age (Fig. 3), with the TPCF successively flattening as the age of the population increases, occupying a larger extent of the disc, and eventually converging to the residual correlation from the large exponential disc distribution in the galaxy.

(iii) We find a range of values of l corr across the sample, from ∼ 100 pc to scales beyond ∼ 2.5 kpc, the largest we can reliably measure given the size of the LEGUS footprint. Similarly, we find a range of fractal dimensions (D 2 ) for young clusters from ∼0.5 to 1.9 across our sample of galaxies (see Table 3). The range of these parameters is substantially larger than the uncertainties, and suggests that there are significant variations in the hierarchical structuring of star clusters from one galaxy to another. Earlier studies show that this is not the case for the gas distribution (see Table 1; Shadmehri & Elmegreen 2011), suggesting that there might be additional physical mechanisms at play in explaining these differences.

(iv) We find signs of some positive correlation of l corr with stellar mass M * , UV-derived star formation rate SFR UV and star formation rate surface density SFR (Fig. 5). We also find relatively stronger and statistically more significant correlation of l corr with the galaxyaveraged Toomre length l T in the disc (Fig. 6), suggesting that rotation-supported gravitational instability might be an important mechanism in setting the scales where gas is hierarchically structured. We stress, however, that we are limited to 12 galaxies in this study, and hence, cannot make fully conclusive inferences.

(v) We demonstrate that we can robustly infer an estimate for the radial scale length of the star cluster distribution in the galaxy (r c ) from the TPCF of its more randomly distributed older clusters, but only for galaxies where the field of view within which we measure star cluster positions is substantially larger than the radial scale length (Fig. 8).

Overall, our results suggest that the hierarchical structure of star clusters, both old and young, is not universal, but instead depends on the physical properties of the host galaxy. For older clusters, this dependence is relatively trivial, since as the cluster population ages, it loses the hierarchical structure with which it formed, and the resulting TPCF simply reflects the overall size of the galaxy. More intriguingly, though, even for young clusters, we measure statistically significant variations in both the fractal dimension and the largest scale of the hierarchical distribution, and show that these correlate with largescale galactic properties. Therefore, cluster formation is possibly not a universal process that operates the same way in all galaxies, which suggests significant scope for future work by extending our study to a larger sample, within which the correlations between cluster distributions and galactic properties of which we see hints can be more reliably measured.

We fit this analytical form to our measured TPCFs and investigate at what point our fitted values of α and A differ from the values for the input, non-truncated fractal distribution by more than 10 per cent. We show our computed TPCFs from the truncated data along with the true TPCFs for each value of R max in Fig. A1. We plot ω instead of 1 + ω in order to make the edge effects more clearly visible.

In general, we find that our TPCFs for the truncated data match analytic expectations to better than 10 per cent for separations x ∼ < R max /5, but that for the truncated-data TPCFs ω falls off much more shallowly than predicted by the analytical relation. This disagreement is likely due to the data-random cross correlation term of the Landy & Szalay (1993) estimator, which becomes dominant at separation close to the size of the field-of-view. Given this result, we set l edge = R max /5, and discard our measured TPCFs at larger separations. However, we caution that our choice R max /5 is somewhat arbitrary, since the divergence between the measured and true TPCFs in our idealized experiment occurs over a finite range of scales, rather than sharply at a single scale.

In Section 4.1, we presented and discussed the TPCF of star clusters divided into young and old clusters based on an age cut (T = 10 Myr). We chose this approach instead of showing the combined TPCF of both young and old clusters, as the physically relevant features in the TPCF are more clearly evident when the sample is divided by age. For completeness, however, we show the combined TPCF in Fig. B1 with their best-fitting models obtained from the methodology outlined in Section 4.2.1 overplotted.

. TPCF of star clusters of all ages for each galaxy in the sample, with their best-fitting models overplotted, using the approach outlined in Section 4.2.1. The best-fitting parameters appropriate to the best-fitting models are denoted on the plot, and the grey shaded regions denote separations where edge effects might play a role, as in Fig. 2.

. Left: The computed TPCF for fractal models with different input 2D fractal dimensions (D 2 ) as indicated in the legend. The shaded region in blue marks separations beyond the maximum scale of the hierarchy (i.e. x > 1/2 l base ). As we can see, the TPCF is a pure power law up to scales where the hierarchy extends, beyond which it sharply flattens, and the slope of the power law is progressively shallower at larger D 2 , with a completely flat TPCF for a purely random distribution (i.e. D 2 = 2.0). Right: Comparison of input fractal dimension provided to the toy model (D 2 ) and the derived 2D fractal dimension α + 2, obtained from the a least-squares linear fit to the TPCF shown in the left-hand panel at x < 1/2 l base . The one-to-one relation is shown as a dashed green line, and error bars indicate the 1σ uncertainties returned by the fit. This shows that the TPCF can reasonably reproduce D 2 from the slope α of the power law.

In this section, we describe a set of physically motivated toy models that we use to infer the features seen in the star cluster TPCFs of the galaxies in our sample in Section 4.1. These three models are meant to characterize the three fitting functional forms described in Section 4.2.1, namely a single power law (Model S), a piecewise power law (Model PW), and a power law with an exponential fall-off (Model PF). We explain the three classes of features by a pure fractal distribution, fractal distribution that transitions to a random one beyond some outer scale, and a radially exponential disc distribution, respectively. Below we discuss the parameters of the model, and how the TPCFs of the model depend on the parameters.

Our procedure for constructing fractal distributions of points uses the same reverse box-counting method previously employed by a number of authors (Bate et al. 1998;Cartwright & Whitworth 2004;Gouliermis et al. 2014;Elmegreen 2018). The method is as follows: we begin with a square of side length L box , which we divide up into 2 l square cells of side length L box /2 l ; where l ≥ 0 is the level in the hierarchy. We start at a base level l base by marking all cells at that level as 'active'. We then subdivide each active cell into four subcells at level l = l base + 1, and randomly decide whether to mark those subcells as active, with probability p = 2 D 2 -1 . We then repeat this procedure recursively: for each active cell at level l, we subdivide it into four cells and level l + 1, which we mark active or inactive with probability p, and so forth. The algorithm terminates at some predetermined maximum level l max ; we place a point in each active cell on this level, with the location of the point set equal to the location of the cell center plus a small random dither to avoid an overly gridded structure.

As described, this algorithm is fully characterized by the following parameters:

(i) D 2 : The 2D fractal dimension of the distribution.

(ii) L box : The maximum spatial extent of the box in which the fractal is present.

(iii) l base : The minimum level, which consequently sets the maximum separation scale L max up to which there is fractal structure.

(iv) l max : The maximum level, which sets the minimum separation L min of the hierarchically (fractal) distributed points.

Thus, in such a setup, we should expect scale-free behaviour in the ranges of separations from L min ∼ L box /2 lmax to L max ∼ L box /2 lmin . For the first set of fractal models we prepare, we vary the input fractal dimension D 2 , and use fixed values of l base = 2 and l max = 14, L box = 1.0 for convenience. The TPCFs for the various input fractal dimensions in such a case are shown in Fig. C1. The TPCF is clearly a pure power law up to L max , beyond which it sharply flattens, and the slope of the power law is progressively shallower at larger D 2 , with a completely flat TPCF for a purely random distribution (i.e. D 2 = 2.0). In addition, we verify that the slope obtained from the power law part of the TPCF for the fractals matches the theoretical prediction, i.e. D 2 = 2 + α (Calzetti et al. 1988;Gomez et al. 1993;Larson 1995). This is shown in the right-hand panel of Fig. C1, where we compare our input fractal dimension to the model D 2 with 2 + α, where we determine α by performing a least-squares fit to the data shown in the left panel at x < 1/8. As we can see, the slopes extracted from the power-law TPCFs match the analytical prediction reasonably well.

In addition to this, we also attempt to vary the maximum scale of the hierarchy L max in our fractal models, keeping the fractal dimension D 2 fixed. This might be important in setting the separation where the TPCF flattens to a value of 1 + ω ∼ 1. This becomes relevant for galaxies whose TPCF is best fit by Model PW, where the slope of the TPCF changes from a steep one to a relatively much shallower one, characteristic of random distributions. We attempt four different values of L max = 1/8, 1/16, 1/32, and 1/64, keeping D 2 = 1.0, L min = 1/2 14 , and L box = 1.0 fixed. We show the results in Fig. C2. We find that the TPCF flattens more or less at the scale of L max , with a fit to the functional form of Model PW yielding a transition point β ≈ L max . This suggests that the transition identified in the power-law behaviour of Model PW is capturing a physical transition in the underlying distribution from one that is fractal/hierarchical in nature, to a mostly random distribution.

Here, we describe the toy models we use to represent the large-scale distribution of star clusters in thin, radially exponential disc inclined at an arbitrary angle relative to the line of sight. Our model contains five parameters:

(i) r c : The exponential scale radius of the distribution (ii) z h : The Gaussian scale height of the galaxy (iii) i: Line-of-sight inclination angle of the galaxy (iv) R max : The maximum radial extent of the galaxy up to which the points are distributed (v) r min : The minimum radius at which points can be found from the center of the galaxy Given these parameters, the model probability density is P (r) = 1 rc exp -r rc ∀ r min ≤ r ≤ R max , P (z) = N (0, z h ), P (θ) = U(0, 2π ), (C1)

where r, z, θ are the coordinates of a cylindrical coordinate system with its origin at the galaxy center and the galaxy mid-plane lying at z = 0, N (μ, σ ) is the Gaussian distribution with mean μ and standard σ , U(a, is a uniform distribution the range b). We generate our galaxy model by drawing (r, z, θ ) coordinates from this distribution, rotating the positions of the points by the chosen inclination angle i, and then de-projecting to obtain the plane-of-sky distribution exactly as we do for observed star clusters (see Section 3.1). We show the resulting TPCFs for a range of values of r c in Fig. C3; we do not show results for varying z h , because we find that the value of this parameter is negligible as long as z h r c . The general shape of the TPCFs is a shallow power law at separations x r c , followed by an fall-off as separations approach x ∼ r c , which is why an exponentially truncated disc is our prime candidate to describe the Model PF fits we obtained in Section 4.2.1. We also attempt to fit a Model PF functional form to the TPCFs of our toy models, and find that it fits very well, with the bestfitting value θ c , i.e. the fitted exponential scale of the fall-off in the TPCF, reproducing the underlying r c quite well. We demonstrate this in the right-hand panel of Fig. C3. In addition, we also attempted introducing an azimuthal pattern, such as a logarithmic spiral, to our exponential disc models. However, we found that the TPCF is relatively insensitive to the introduction of spiral arms, as we show in Fig. C3, apart from a slight excess at the smallest separations.

Figure Left: TPCFs with the axi-symmetric thin exponential disc toy models described in Section C2 for values of r c = 0.05, 0.1, 0.2, and 0.4. The other model parameters are kept fixed at z h , i, R max , r min = 0.2, 0.05, 30 • , 1.0, 0.01. We also show the TPCF for an exponential disc with r c = 0.2 containing logarithmic spiral arms (black dashed), and find that it is more or less identical to that of an axi-symmetric disc, apart from a slight excess of correlation at the smallest separations. Right: The value of θ c we obtain by performing a least-squares fit of the functional form for Model PF (equation 9) to the measured TPCFs for a range of exponential disc scale radii r c . A one-to-one relation is plotted to guide the eye. We find that the parameter θ c reproduces the underlying r c of the distribution quite well.

Comparison of the TPCFs (solid lines with error bars) and fits to θ c (dashed vertical lines) of galaxy disc models with r c , z h , i, r min = 0.2, 0.05, 30 • , 0.05, and three different values of R max = 0.4, 0.6, and 1.0, corresponding to 2r c , 3r c , and 5r c , respectively. As we can see, the fitted value of θ c for the latter case is reasonably close to the true value of 0.2, whereas θ c (∼0.6) for the former two cases overestimates the true value by a factor ∼3. Thus, insufficient radial sampling of the galaxy can lead to an overestimated value for the scale length inferred from the TPCF using Model PF.

The TPCFs for this model also depend weakly on the other parameters of the toy model. We will not discuss these variations further, except to note that the dependence on R max becomes relevant to the discussion in Section 4.4.3, where we find that the inferred scale radii from the Model PF fit to θ c for the larger spiral galaxies is overestimated by a factor 2-3 as compared to other estimates in the literature. We understand this to arise due to the fact that a smaller extent of the entire radial distribution of the clusters would be sampled by the HST field of view for larger galaxies. To test whether this limited field of view can lead to an overestimate of θ c > r c when fitting Model PF, we set up an exponential disc distribution with r c , z h , i, r min = 0.2, 0.05, 30 • , 0.05 and three different values of R max = 0.4, 0.6, and 1.0, which corresponds to 2r c , 3r c , and 5r c , respectively. We then calculate their TPCFs and compare the value of θ c with the input r c . This analysis is shown in Fig. C4. As we can see, when R max = 5r c , the fitted θ c is reasonably close to the input r c = 0.2, whereas for R max = 2r c and 3r c , the fitted θ c is considerably higher by a factor ∼3. This shows that if a galaxy with a given scale length r c is not observed to sufficiently large radii, r ∼ 5r c , then the estimate for θ c will overestimate the true value of r c by factors of a few.

Here we list the values and sources for the physical quantities we used, namely the galaxy-averaged gas surface density g and flat rotational velocity v rot , in the calculation of the average Toomre length in a galaxy l toomre using equation ( 12). For g , we use surface densities of molecular gas as it is the phase of the ISM where star formation is expected to occur (Bigiel et al. 2008). For v rot , we use H I rotation curves available in the literature, as this is the most widely available line that traces the rotational velocities in spiral galaxies. We list the values and references for the galaxies below. Note that in some cases where the source does not report a value of g , we explicitly calculate g by averaging the total molecular gas mass M H2 reported in the source in a disc of radius R 25 , using the values of R 25 given in Table 1. In addition, l toomre is only computed for the spiral galaxies in our sample.

g : Direct estimate: NGC 0628, NGC 5194, NGC 5457, and NGC 6503 from Kennicutt (1998). Indirect calculation: M H2 of NGC 3344 NGC 3627 from Young et al. (1989), M H2 of NGC 1566 from Bajaja et al. (1995), and M H2 of NGC 7793 from Israel, Tacconi & Baas (1995).

v rot : NGC 0628, NGC 3627, and NGC 5194 (THINGS survey, de Blok et al. 2008), NGC 1313 and NGC 7793 (Local Volume H I survey, Wang et al. 2017;Koribalski et al. 2018), NGC 1566 (WAL-LABY, Elagali et al. 2019), NGC 3344 (Meidt et al. 2009), NGC 5457 (Guélin & Weliachew 1970), and NGC 6503 (Greisen et al. 2009).

This has been typeset from T E X/L A T E X file prepared by the author.