# The VMC survey - XLVII. Turbulence-Controlled Hierarchical Star Formation in
the Large Magellanic Cloud

Amy E. Miller,1,2 Maria-Rosa L. Cioni,1 Richard de Grijs,3,4 Ning-Chen Sun,5
Cameron P. M. Bell,1 Samyaday Choudhury,3,4 Valentin D. Ivanov,6 Marcella
Marconi,7 Joana Oliveira,8 Monika Petr–Gotzens,6 Vincenzo Ripepi,7 and Jacco
Th. van Loon8
1Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, D-14482
Potsdam, Germany
2Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-
Liebknecht-Str. 24/25, D-14476 Golm (Potsdam), Germany
3School of Mathematical and Physical Sciences, Macquarie University, Balaclava
Road, Sydney, NSW 2109, Australia
4Research Centre for Astronomy, Astrophysics and Astrophotonics, Macquarie
University, Balaclava Road, Sydney, NSW 2109, Australia
5Department of Physics and Astronomy, University of Sheffield, Hicks Building,
Hounsfield Road, Sheffield S3 7RH, UK
6European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching
bei München, Germany
7INAF – Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, I–80131
Naples, Italy
8Lennard-Jones Laboratories, Keele University, ST5 5BG, UK
E-mail<EMAIL_ADDRESS>(AEM)

(Accepted ?, Received ?; in original form ?)

###### Abstract

We perform a statistical clustering analysis of upper main-sequence stars in
the Large Magellanic Cloud (LMC) using data from the Visible and Infrared
Survey Telescope for Astronomy survey of the Magellanic Clouds. We map over
2500 young stellar structures at 15 significance levels across $\sim$120
square degrees centred on the LMC. The structures have sizes ranging from a
few parsecs to over 1 kpc. We find that the young structures follow power-law
size and mass distributions. From the perimeter–area relation, we derive a
perimeter–area dimension of 1.44$\pm$0.20. From the mass–size relation and the
size distribution, we derive two-dimensional fractal dimensions of
1.50$\pm$0.10 and 1.61$\pm$0.20, respectively. We find that the surface
density distribution is well-represented by a lognormal distribution. We apply
the Larson relation to estimate the velocity dispersions and crossing times of
these structures. Our results indicate that the fractal nature of the young
stellar structures has been inherited from the gas clouds from which they form
and that this architecture is generated by supersonic turbulence. Our results
also suggest that star formation in the LMC is scale-free from 10 pc to 700
pc.

###### keywords:

stars: early type – methods: statistical – stars: formation – galaxies:
individual – galaxies: stellar content – galaxies: individual: Large
Magellanic Cloud – galaxies: structure

††pagerange: The VMC survey - XLVII. Turbulence-Controlled Hierarchical Star
Formation in the Large Magellanic Cloud–B††pubyear: 2021

##  1 Introduction

Stars are observed in galaxies both as individual field stars and grouped in
clusters, associations and complexes. Young stars with ages less than 100 Myr
are more spatially clustered than evolved, older stars. These spatially
clustered young star groups are hierarchically organised: larger, lower-
density groups contain one or more smaller, higher-density groups, which then
fragment into even smaller and more compact groups (Elmegreen, 2010). This
hierarchical architecture is also observed for gas structures in the
interstellar medium (ISM; Elmegreen & Falgarone, 1996; Blitz & Williams,
1999). Observations and simulations suggest that the organisation of young
stellar structures originates from the fractal distribution of the ISM
(Elmegreen & Falgarone, 1996; Efremov & Elmegreen, 1998; Elmegreen &
Elmegreen, 2001; Vázquez-Semadeni et al., 2016). This is evidence of
hierarchical star formation, also known as ‘scale-free’ or fractal star
formation (Gomez et al., 1993; Larson, 1995; Bate et al., 1998).

Classically, it was thought that stars could form in either the field of a
galaxy as single stars or concurrently in star clusters that formed from the
same molecular cloud. This would mean that there are two distinct star
formation modes: ‘clustered’ and ‘field’ or ‘distributed’ star formation (Lada
& Lada, 2003). However, Bressert et al. (2010) found that young stellar
systems exist in a continuous range of stellar densities or masses, implying
that there is no bimodality in star formation and star clusters only represent
the densest part of the continuum. Further investigation of the hierarchical
properties of star formation is necessary to gain a better understanding of
the full continuum over a large parameter range and assess whether there are
any characteristic scales or distinct modes of star formation, which may
correspond to some physical mechanisms.

Hierarchical star formation has many broad-reaching astrophysical
implications. For example, the universal surface brightness profiles of young
massive clusters in different environments (Grudić et al., 2018) can be
explained if clusters form through the hierarchical merging of subclusters. In
addition, hierarchical mergers of subclusters may lead to a phase when the
cluster becomes very compact, characterised by short dynamical time-scales,
and dynamical mass segregation may occur during this very early stage (Allison
et al., 2009). Binary star clusters that are located very close to each other
and have similar ages (e.g., the Double Cluster; Keller et al., 2001) may have
formed as a consequence of hierarchical star formation and may be affected by
mutual dynamical interactions. Additionally, the commonly observed multiple
stellar populations in star clusters could be a consequence of hierarchical
star formation: a newly formed star cluster may continue to accrete gas from
its surroundings and form new stars of younger ages (Bekki, 2017). On a larger
scale, the observed ultra-compact dwarf galaxies may be the result of mergers
of many star clusters in a star-forming complex (Urrutia Zapata et al., 2019),
which is expected in hierarchical star formation. Hierarchical star formation
can also explain why only a minority of core-collapse supernovae are found in
star clusters (e.g. Sun et al., 2020) while the majority are not.

The properties of the star-forming hierarchy are not yet fully understood
since they have only been studied for a limited number of star-forming regions
and galaxies. The goal of this paper is to better understand the star-
formation hierarchy on scales between star clusters and galaxies. Performing a
global study of hierarchical star formation in the Milky Way is very
difficult, because we cannot view the Galaxy in its entirety and there is dust
obscuration in the Galactic plane. This, therefore, turns our attention to the
Milky Way’s most massive satellite galaxy, the Large Magellanic Cloud (LMC).
The LMC is seen face-on and features active star formation. It is the ideal
candidate to study hierarchical star formation owing to its orientation,
proximity and the small line-of-sight depth (ranging from 2.6 to 4.2 kpc in
the disc and Bar regions; Subramanian & Subramaniam, 2009). The LMC exhibits
signatures of interactions with its companion, the Small Magellanic Cloud
(SMC), as well as with the gravitational potential and halo gas of the Milky
Way. The LMC was the first nearby galaxy in which a large population of young,
massive stars was observed. Hodge (1973) observed star clusters in the LMC and
noted that their distribution is connected to the patchy distribution of Hii
regions. They also noted a pattern of star formation on 1 kpc size scales,
which is related to very dense gas regions in the LMC. Efremov & Elmegreen
(1998) showed that the average age difference between pairs of star clusters
in the LMC increases with their separation. This suggested that star formation
is hierarchical in space and time. They showed that, in small regions, stars
form quickly, whereas in large regions, which also contain smaller groups of
stars, stars form over a longer period. Bastian et al. (2009) analysed the age
distributions of different stellar groups in the LMC. They found groups of
younger stars to be highly substructured and nonuniform, and that this
nonuniform distribution dissolves into the background, uniform stellar
distribution on a time-scale of $\sim$175 Myr. Bonatto & Bica (2010) found
that the young star clusters from the catalogue of Bica et al. (2008) are
spatially correlated with one another and with star-forming structures in the
ISM. However, the catalogue of Bica et al. (2008) was created from earlier
catalogues originating from different studies, thus making selection biases
and incompleteness difficult to judge.

Using data from the Visible and Infrared Survey Telescope for Astronomy
(VISTA) Magellanic Clouds survey (VMC; Cioni et al., 2011), Sun et al. (2017a,
b, 2018) found that young stellar structures in the Magellanic Clouds have
power-law size and mass distributions which can be best explained by a
scenario of hierarchical star formation. Zivkov et al. (2018) identified 31
pre-main-sequence (PMS) structures in a 1.5 deg2 region just west of the
Tarantula Nebula and north of the Bar in the LMC, based on VMC data, with a
mass distribution slope which was consistent with hierarchical star formation
regulated by turbulence.

Sun et al. (2017a, b, 2018) used a contour-based clustering technique to
identify the projected two-dimensional (2D) boundaries of groupings of young
stars. The technique was first used in nearby galaxies (Gouliermis et al.,
2010; Gouliermis et al., 2015, 2017). This method works by identifying stellar
overdensities in a sequence of significance levels above the background
stellar density from a surface density map of young stars. It is particularly
efficient at partitioning spatial hierarchies. This technique has never been
applied to the entire ($\sim$120 square degrees) LMC.

To better understand the star-formation hierarchy, in this paper we expand the
work begun by Sun et al. (2017a, b) and exploit the VMC survey’s deep near-
infrared photometry and high spatial resolution ($<1\arcsec$) to perform a
galaxy-wide, contour-based clustering analysis of young stellar structures
across the LMC. The VMC photometry in the near-infrared $YJK_{\text{s}}$ bands
is less affected by extinction effects than photometry in optical passbands.
Following Sun et al. (2017a, b, 2018), we will refer to all stellar groups,
clusters, associations and complexes as young stellar structures.

In this paper we identify young stellar structures on size scales from a few
parsecs to over 1 kpc (bridging from cluster scales to galaxy scales) in the
LMC, estimate their parameters and compare with past studies of hierarchical
star formation to understand the physical mechanisms behind this seemingly
universal architecture. This paper is organised as follows. In Section 2, we
describe the VMC data used, the sample selection of young stars and our
identification of statistically grouped young stellar structures. In Section
3, we show a dendrogram of the second-largest structure identified, we
calculate parameters of all young structures, and analyse relevant
correlations and distributions. We compare with past studies in the Magellanic
Clouds, nearby galaxies and the ISM in Section 4. Section 5 contains our
conclusions.

##  2 Data and Young Stellar Structure Identification

###  2.1 VMC Data

Figure 1: The 68 VMC LMC tiles we study in this paper, overlaid on a stellar
density map with sources characterised by $J$ and $K_{\rm s}<20$ mag. The
equatorial coordinates were transformed using a zenithal equidistant
projection with the origin at ($80^{\circ}.98,-69^{\circ}.69$), which
corresponds to the dynamical centre of young carbon stars in the LMC (Wan et
al., 2020).

The VMC survey (Cioni et al., 2011) is a near-infrared $YJK_{\text{s}}$ survey
of the Magellanic System covering 170 deg2. The data were gathered between
2009 and 2018 with the 4.1 m VISTA telescope (Sutherland et al., 2015),
operated by the European Southern Observatory (ESO), using the VISTA infrared
camera (VIRCAM; Dalton et al., 2006). VIRCAM has 16 detectors arranged in a
$4\times 4$ array; one exposure is a ‘pawprint’. To cover the gaps between the
detectors and ensure contiguous coverage, six pawprints shifted by a specific
offset were combined to form a VMC ‘tile’. Each VMC tile covers $\sim$1.77
deg2. The pawprints are available at both the VISTA Science
Archive111http://horus.roe.ac.uk/vsa/ (VSA; Cross et al., 2012) and the ESO
Science Archive Facility222http://archive.eso.org/cms.html. They were reduced
using the VISTA Data Flow System (Irwin et al., 2004; González-Fernández et
al., 2018), version 1.5, by the Cambridge Astronomy Survey Unit (CASU). Point-
spread-function (PSF) photometry was carried out on homogenised, stacked tile
images and photometric errors and local completeness have been calculated
based on artificial star tests (for more details about the PSF photometry, see
Rubele et al., 2012, 2015). We downloaded the VMC photometric catalogue data
from the VSA. Although the VISTA images are publicly available, the PSF
photometry used in our work is proprietary to the VMC team, with a data
release planned in early 2022.

There are 68 VMC tiles (shown in Figure 1) covering the LMC, and we include
them all in this paper. There are a total of 116,336,429 sources in the 68
tiles. First, we selected stars with $J$ and $K_{\rm s}<22$ mag and
photometric errors $<0.15$ mag in either filter. After this selection step
there were 38,802,981 sources. We next cross-matched adjacent tiles and
removed duplicate sources located within 0.5″of one another, leaving us with
36,030,318 sources. We make additional selection cuts in the next section.
Figure 1 shows a vertical gap in the stellar distribution. This gap is owing
to an observational mistake in which pawprints of tile LMC 4$\\_$4 were
shifted too far eastwards (in Right Ascension). Complementary VISTA
observations are planned with the goal of filling the gap. In this paper, we
fill the vertical gap using optical data from the Survey of MAgellanic Stellar
History (SMASH; Nidever et al., 2017, 2021). This process is described in
Appendix A.

####  2.1.1 Selection of VMC Upper Main-Sequence Stars

Figure 2: CMDs of VMC tiles LMC 4$\\_$6 (top row) and LMC 6$\\_$6 (bottom
row). All isochrones plotted are of ages $\log(\tau\text{ yr}^{-1})=$ 7.0, 8.0
and 9.0, in black, from left to right. (Top left) CMD of tile LMC 4$\\_$6 with
magnitudes uncorrected for extinction, and with isochrones and the UMS
selection box overlaid in red. (Top middle) CMD of tile LMC 4$\\_$6 with
dereddened magnitudes, isochrones and the UMS selection box overlaid in red.
(Top right) CMD of tile LMC 4$\\_$6 with dereddened magnitudes and components
overlaid. (Bottom) As the top row, but for tile LMC 6$\\_$6\.

Figure 2 shows $J-K_{\rm s}$ versus $K_{\rm s}$ colour–magnitude diagrams
(CMDs) of VMC tiles LMC 4$\\_$6 (top panels) and 6$\\_$6 (bottom panels). We
use this colour combination because it is less affected by reddening than
($Y-K_{\rm s}$). We estimate the prevailing extinction based on the Magellanic
Clouds’ reddening map provided by Skowron et al. (2021), derived from the
properties of red clump (RC) stars and presented in units of $E(V-I)$, with a
resolution of $1\arcmin.7\times 1\arcmin.7$ in the central regions of the LMC.
We used $A_{I}=1.5E(V-I)$ (Skowron et al., 2021) and the calibrations of
Schlegel et al. (1998)’s $E(B-V)$ values – which we refer to as
$E(B-V)_{\text{SFD}}$ – provided by Schlafly & Finkbeiner (2011), to express
the extinction coefficients in terms of $E(V-I)$. The Schlafly & Finkbeiner
(2011) calibration we used is $A_{I}=1.505E(B-V)_{\text{SFD}}$. Using these
two relations and Appendix B of González-Fernández et al. (2018), we derived
$A_{J}=0.703E(V-I)$ and $A_{K_{\text{s}}}=0.307E(V-I)$. We then dereddened
each stellar source by referring to the closest point in the Skowron et al.
(2021) map.

In Figure 2 we also overlaid PARSEC version 1.2S (Bressan et al., 2012)
isochrones of ages $\log(\tau\text{ yr}^{-1})=$ 7.0, 8.0 and 9.0 and a
metallicity of [Fe/H] = $-0.3$ dex, which is typical of young OB stars in the
LMC (Rolleston et al., 2002). In each tile we analyse, the isochrones are
shifted by the LMC’s distance modulus, $(m-M)_{0}=18.49$ mag (de Grijs et al.,
2014). We did not account for the inclination of the LMC’s disc and used the
same distance modulus for all tiles; accounting for the inclination of the
disc would slightly shift the position of the isochrones vertically. The
isochrones were also shifted from the Vega system to the VISTA system using
the relations of González-Fernández et al. (2018):
$J_{\text{VISTA}}=J_{\text{Vega}}$ and
$K_{\text{s}_{\text{VISTA}}}=K_{{\text{s}}_{\text{Vega}}}-0.011$.

The third panel in each row of Figure 2 shows the different components visible
in the near-infrared CMDs, including Milky Way foreground stars, the main
sequence (MS), the RC, the red-giant branch (RGB), red supergiants (RSGs) and
background galaxies. To select young MS stars, we looked at the bright young
MS, and therefore we selected upper-MS (UMS) stars. We adopted selection
criteria based on the PARSEC evolutionary models and similar to those of Sun
et al. (2017a, b). We selected UMS stars if they fell into the colour range
$-0.3<J-K_{\rm s}<0.19$ mag. This was done to differentiate UMS stars from RGB
and RC stars. We chose this wider colour range similarly to Sun et al.
(2017b), because although we correct for extinction, the ISM is inhomogeneous
across the LMC, which leads to an increase in the width of the MS. We also
adopted a magnitude range of $13.0<K_{\rm s}<18.1$ mag. The brighter limit at
13 mag allows for some blue-loop stars to fall into our selection criteria,
but stars do not spend much time at this evolutionary stage and contamination
is therefore minimal. The fainter limit at 18.1 mag was chosen to ensure that
all of our UMS stars are younger than 1 Gyr. This is evidenced by the oldest
isochrone of 1 Gyr; see Figure 2. After making these cuts, 362,047 sources
remain.

Next, we selected stars with photometric errors less than 0.10 mag, which left
us with 355,250 sources. Finally, we made cuts based on the stellar
probability parameter, $\tt{Star\\_prob}$, and the $J$ and $K_{\text{s}}$
sharpness indices, $\tt{SHARP}$. As Bell et al. (2019) showed (their figure
2), stellar objects should have $\tt{Star\\_prob}\geq 0.34$ and $J$ and
$K_{\text{s}}$ sharpness indices $\tt{SHARP}\leq 0.5$. After making these
cuts, our selection retained 307,581 UMS stars.

In the region of our selection box shown in Figure 2, contamination owing to
foreground Milky Way stars is very small, since they are generally found at
colours around $(J-K_{\text{s}})\sim 0.7$ mag. We checked the contamination by
Galactic sources by cross-matching the VMC LMC sources with Gaia early Data
Release 3 (Gaia Collaboration et al., 2021). We found in our UMS region that
the parallaxes for all matches were $\ll 0.2$ mas, which indicates that
contamination by the Milky Way is negligible in our colour–magnitude region
(for a more in-depth discussion of foreground contamination in the VMC data,
see El Youssoufi et al., 2021). Additionally, we checked for contamination
using the TRIdimensional modeL of thE GALaxy (TRILEGAL) code (Girardi et al.,
2005). This code estimates the number of Milky Way stars within a projected
sky area. We ran the code on an area of 10 square degrees centred on the LMC.
The simulation showed that there is 0.2% Milky Way contamination in our
selection box. Therefore, we conclude that Milky Way contamination is indeed
negligible.

###  2.2 Stellar Surface Density Map

Figure 3: (Left) Surface density map of UMS stars selected from the VMC data.
The colour bar is given in units of stars pc-2. Labeled are a number of LMC
features visible in the young stellar population, as well as three star-
forming regions. (Right) Detected young stellar structures coloured by their
significance level. The right-hand colour bar is labelled by the significance
level and the number of structures contained in each level. The two largest
structures are labelled $S_{1}$ and $S_{2}$.

To explore the spatial distribution of UMS stars, we applied kernel density
estimation (KDE) – convolution of the stellar distribution with a Gaussian
kernel – to a 2D binned distribution of the UMS sample. In our application of
KDE, the kernel size specifies the resolution of the map. The choice of kernel
size is arbitrary and depends on the image resolution, the distance to the
galaxy and the size of the objects considered. It is best to experiment to
find the optimal kernel size (Gouliermis et al., 2017; Sun et al., 2018). The
resolution of the VMC survey is 1 arcsec, which corresponds to 0.24 pc at the
distance of the LMC. We first grid our UMS into bins of 1 pc $\times$ 1 pc,
chosen to be slightly larger than the VMC resolution limit. Next, we
experimented with kernel sizes from 5 pc to 20 pc. We found a minimum possible
kernel of 10 pc. Sun et al. (2017b, 2018) also found that a kernel of 10 pc
represents an optimal balance between resolution and noise. Using these
constraints, we created a stellar surface density map. Based on our
resolution, we can confidently identify structures larger than 10 pc. The left
panel of Figure 3 shows the KDE map with three prominent components of the LMC
seen in the young stellar population: the Bar, the Northwestern Arm and the
Southwestern Arm. Also pictured are three well-known star-forming regions: 30
Dor, Shapley Constellation III (S.III) and N11; see El Youssoufi et al. (2019)
for a more in-depth discussion of the morphological features of the LMC’s
young population. We filled the gap near tile LMC 4$\\_$4 with optical SMASH
data; see Appendix A. The mean surface density of the KDE map is 0.003 stars
pc-2, the median surface density is $4\times 10^{-5}$ stars pc-2 and the
standard deviation is 0.006 stars pc-2. These values are similar to those
found by Sun et al. (2018) for the SMC. Most of the values in our map (see the
left-hand panel of Figure 3) and in Sun et al. (2018)’s map are near zero, but
there are some outliers which cause the mean density to be higher than the
median.

###  2.3 Young Stellar Structure Detection

Young stellar structures are identified from the KDE map (see the left-hand
panel of Figure 3) for a range of significance levels above the median density
of the KDE map: from 1$\sigma+$ the median stellar KDE density to 15$\sigma+$
the median stellar KDE density, in steps of 1$\sigma$. The median stellar
density of the map is effectively zero; therefore, using the median does not
affect the statistics of the less populous regions of the LMC, whereas using
the mean would. This process resulted in the detection of more than 7000
structures at the 15 significance levels we adopted. To remove spurious
detections, we applied two selection criteria: (1) each candidate structure
must enclose at least five stars, and (2) every candidate structure on the
1$\sigma$ and 2$\sigma$ levels must contain one or more contours at higher
significance levels. The limit of five stars was chosen to coincide with that
used by Sun et al. (2017b, 2018), and it has also been used in similar earlier
studies to reduce contamination (Bastian et al., 2007; Bastian et al., 2009;
Gouliermis et al., 2017). After application of these criteria, we identified a
total of 2846 structures at the 15 significance levels. The result is shown in
the right-hand panel of Figure 3; the legend on the right shows the different
significance levels (1$\sigma$ to 15$\sigma$) and below each level the number
of structures identified at each level is indicated. In Appendix B we compare
our catalogue of young stellar structures with similar catalogues from the
literature and explain any differences we find. We also estimate that 45% of
our structures were not previously found.

What is the definition of a stellar structure or stellar grouping in the
context of our study? As we show in the next sections, our identified
structures range in radius from a few parsecs to over 1 kpc. Therefore, our
structures include individual stars, young clusters, associations and
complexes. The largest identified structures might be associated with galaxy-
scale structures or processes. We note that it is difficult to differentiate
among clusters, associations, complexes and larger stellar groups. We cannot
assess if a group of identified young stars is gravitationally bound, nor
whether the stellar group members are located truly close in space. Therefore,
if our algorithm identifies a group above the median surface density which
passes both of our selection criteria, we refer to it as a young stellar
structure or group.

The right-hand panel of Figure 3 displays the nonuniform design and
organisation of the detected young stellar structures. The hierarchical
organisation is apparent, especially considering the two largest structures,
which we label $S_{1}$ and $S_{2}$. To further demonstrate the hierarchical
structure, we created a dendrogram of $S_{2}$, shown in Figure 4. The
dendrogram shows the parent–child links between the young stellar structures
at the 15 significance levels.

Figure 4: Dendrogram of $S_{2}$, the second-largest structure we identified,
showing the corresponding parent–child relationships. $S_{2}$ has 451
structures contained within the 14 higher significance levels.

##  3 Results

###  3.1 Parameters of the Young Stellar Structures

In the following subsections, we show how we derived structural parameters
using the boundaries of the structures shown in the bottom panel of Figure 3.
The basic parameters are the structures’ radii $R$ (in pc), perimeters $P$ (in
pc) and areas $A$ (in pc2). We account for stellar completeness and calculate
the completeness-corrected numbers of stars in each structure, $N_{*}$. We
also calculate the surface density, $\Sigma$, of each structure. Finally, we
use the Larson (1981) relation to estimate velocity dispersions ($\sigma_{v}$)
and then use it to derive different crossing times ($t_{\text{cr}}$).

####  3.1.1 Fundamental Parameters

All parameters of the young stellar structures are based on the 2D projected
boundaries of the structures identified in Section 2.3. Based on the vertices
of each structure, we calculated their areas and perimeters. The radius of
each structure is defined as the radius of a circle with the same area as that
delimited by each contour.

In order to accurately estimate the number of stars in each structure, we
corrected for photometric incompleteness. We followed the same procedure as
Sun et al. (2017b): we first assigned weights to each star,
$w=1.0/\text{min}[f_{J},f_{K_{\rm s}}]$, where $f_{J}$ and $f_{K_{\rm s}}$ are
the local completeness levels in the $J$ and $K_{\rm s}$ bands, respectively.
The completeness parameters were derived by adding artificial stars to the
images, then searching for them in the photometric catalogues (see Rubele et
al., 2012, 2015). Next, we calculated the completeness-corrected number of
stars, $N_{*}$, with

$N_{*}=\sum_{i}w_{i}-\Sigma_{\text{bg}}A,$ (1)

where $i$ accounts for each star in the contour, $\Sigma_{\text{bg}}$ is the
median density of the KDE map and $A$ is the area within the contour. The
first term corrects for photometric incompleteness and the second term
estimates the number of background stars in chance alignment.

It has been shown by Sun et al. (2018) that the number of selected UMS stars
in young stellar structures is proportional to their total masses to a good
approximation (their figure 12). Therefore, we also used $N_{*}$ as a proxy
for mass. Finally, we calculated the stellar surface density,
$\Sigma=N_{*}/(\pi R^{2})$.

####  3.1.2 Velocity dispersion and crossing time estimated from the Larson
relation

Larson (1981) identified a power-law scaling relation between velocity
dispersion and molecular cloud size:

$\sigma_{v}\approx 1.1(2R)^{0.38}.$ (2)

This relation is known as the Larson relation (for an updated discussion on
the exact form of the Larson relation, see Cen, 2021). It is similar to the
Kolmogorov law for incompressible turbulence and is therefore expected to
result from turbulence. Studies of young stars in star-forming regions in the
Milky Way also show that young groups of stars are generally consistent with
the Larson relation (Ha et al., 2021). We used the Larson relation to estimate
the total three-dimensional (3D) velocity dispersions of the r.m.s. velocity
of the internal stellar motions in each structure, $\sigma_{v}$.

The crossing time is defined as the time it takes a star to travel a distance
equivalent to the size of a structure’s half-mass size (Elmegreen, 2000;
Gouliermis et al., 2017; Ha et al., 2021):

$t_{\text{cr}}\approx\frac{2R}{\sigma_{\text{v}}}.$ (3)

We calculated the crossing times for all identified stellar structures using
Equation 3. The identified structures are listed in Table 1. Included are
their IDs, significance levels, geometric center coordinates, and physical
parameters. Only the first 10 structures are shown, the rest are available
online.

Table 1: The identified young stellar structures and their physical parameters. Columns include ID, significance level ($\sigma$), equatorial coordinates corresponding to the geometric center ($\alpha$(J2000) and $\delta$(J2000) in deg), the radius ($R$ in pc), number of stars ($N_{*}$), surface density ($\Sigma$ in pc-2), velocity dispersion ($\sigma_{v}$ in km s-1), and crossing time ($t_{\text{cr}}$ in Myr). The first 10 entries are shown here as an example; the complete catalogue is available online. ID | Level |  $\alpha$(J2000) |  $\delta$(J2000) | $R$ | $N_{*}$ | $\Sigma$ | $\sigma_{v}$ | $t_{\text{cr}}$
---|---|---|---|---|---|---|---|---
| ($\sigma$) | (deg) | (deg) | (pc) |  | (pc-2) | (km s-1) | (Myr)
1 | 1 | 93.214 | -68.268 | 40.9. | 121.1 | 0.023 | 5.6 | 14.6
2 | 1 | 90.320 | -67.738 | 26.0. | 29.0 | 0.014 | 4.7 | 11.0
3 | 1 | 90.021 | -68.635 | 24.4. | 36.6 | 0.019 | 4.6 | 10.5
4 | 1 | 89.726 | -68.523 | 41.6. | 113.0 | 0.021 | 5.6 | 14.7
5 | 1 | 89.556 | -68.297 | 22.4. | 17.1 | 0.011 | 4.5 | 10.0
6 | 1 | 89.452 | -68.464 | 34.1. | 60.6 | 0.017 | 5.2 | 13.0
7 | 1 | 89.523 | -68.624 | 26.9. | 36.9 | 0.016 | 4.8 | 11.2
8 | 1 | 89.377 | -69.194 | 36.9. | 121.0 | 0.028 | 5.4 | 13.7
9 | 1 | 84.010 | -67.089 | 1008.6 | 33832.9 | 0.011 | 18.4 | 109.8
10 | 1 | 88.954 | -68.205 | 34.9 | 49.9 | 0.013 | 5.3 | 13.2
Table 2: Demographics of the young stellar structures at each significance level. Columns include information about the significance level ($\sigma$), number of structures ($N$), minimum radius ($R_{\text{min}}$), average radius ($\langle R\rangle$), maximum radius ($R_{\text{max}}$), average number of stars ($\langle N_{*}\rangle$), average surface density ($\langle\Sigma\rangle$), average magnitude of the fifth brightest star ($\langle K_{\text{s0}}\rangle$), average velocity dispersion ($\langle\sigma_{v}\rangle$) and average crossing time ($\langle t_{\text{cr}}\rangle$) of the young stellar structures at each level. Level | $N$ | $R_{\text{min}}$ | $\langle R\rangle$ | $R_{\text{max}}$ | $\langle N_{*}\rangle$ | $\langle\Sigma\rangle$ | $\langle\ K_{\text{s0}}\rangle$ | $\langle\sigma_{v}\rangle$ | $\langle t_{\text{cr}}\rangle$
---|---|---|---|---|---|---|---|---|---
($\sigma$) | (structures) | (pc) | (pc) | (pc) | (stars) | (stars pc-2) | (mag) | (km s-1) | (Myr)
1 | 48 | 19 | 120 | 1839 | 4422$\pm 23797$ | 0.01$\pm$0.00 | 16.51 | 6.7$\pm$3.5 | 22.6$\pm$25.9
2 | 233 | 11 | 40 | 1198 | 667$\pm 7509$ | 0.02$\pm$0.01 | 16.93 | 5.1$\pm$1.6 | 13.0$\pm$9.6
3 | 701 | 3 | 19 | 676 | 164$\pm 2115$ | 0.04$\pm$0.02 | 16.76 | 4.0$\pm$1.1 | 8.4 $\pm$5.4
4 | 476 | 3 | 18 | 488 | 170$\pm 1592$ | 0.05$\pm$0.02 | 16.56 | 3.8$\pm$1.2 | 8.0 $\pm$5.4
5 | 317 | 2 | 20 | 232 | 173$\pm 752$ | 0.06$\pm$0.05 | 16.23 | 4.0$\pm$1.2 | 8.5 $\pm$5.2
6 | 273 | 3 | 18 | 162 | 139$\pm 491$ | 0.08$\pm$0.04 | 16.12 | 3.8$\pm$1.1 | 7.9 $\pm$4.5
7 | 196 | 2 | 18 | 106 | 135$\pm 310$ | 0.10$\pm$0.08 | 15.79 | 3.8$\pm$1.1 | 8.0 $\pm$4.2
8 | 161 | 3 | 16 | 86 | 113$\pm 221$ | 0.11$\pm$0.05 | 15.69 | 3.8$\pm$0.9 | 7.7 $\pm$3.4
9 | 118 | 4 | 15 | 78 | 110$\pm 184$ | 0.12$\pm$0.04 | 15.45 | 3.7$\pm$0.9 | 7.5 $\pm$3.1
10 | 90 | 2 | 15 | 70 | 107$\pm 164$ | 0.14$\pm$0.08 | 15.33 | 3.7$\pm$0.9 | 7.3 $\pm$3.0
11 | 68 | 3 | 14 | 44 | 106$\pm 116$ | 0.15$\pm$0.09 | 15.17 | 3.7$\pm$0.8 | 7.2 $\pm$2.6
12 | 54 | 3 | 14 | 32 | 101$\pm 95$ | 0.16$\pm$0.05 | 15.08 | 3.6$\pm$0.7 | 7.0 $\pm$2.4
13 | 44 | 3 | 13 | 27 | 97$\pm 84$ | 0.18$\pm$0.07 | 14.99 | 3.6$\pm$0.7 | 6.8 $\pm$2.2
14 | 38 | 3 | 11 | 26 | 87$\pm 76$ | 0.22$\pm$0.13 | 14.11 | 3.4$\pm$0.7 | 6.4 $\pm$2.0
15 | 29 | 3 | 11 | 25 | 87$\pm 74$ | 0.21$\pm$0.09 | 14.08 | 3.4$\pm$0.6 | 6.4 $\pm$1.7

####  3.1.3 Demographics

Table 2 provides an overview of the average parameters at each level. From the
lowest (1$\sigma$) to the highest levels (15$\sigma$), the structures’ average
parameters change: the average $R$, $N_{*}$ and the average $K_{\text{s}}$
magnitude of the fifth brightest star decrease, whereas the average $\Sigma$
increases. The velocity dispersions are highest and with the largest standard
deviations at the 1,2$\sigma$ levels and lowest and more centred around the
mean at the highest significance levels. The crossing times are longest at the
1,2$\sigma$ levels and have smaller standard deviations at the lowest levels.
These trends demonstrate the hierarchical nature of the structures, implying
that the large lower-density structures contain one or more higher-density
groups.

Estimating the ages of each structure is complicated, and we do not attempt to
derive them in this paper. However, we can compare with similar regions in
VMC-CMD space defined for the Magellanic Clouds by El Youssoufi et al. (2019).
They defined three regions that overlap with our selection box: (1) region A,
a population of stars with a median age of $\sim$20 Myr; (2) region B, a
population of stars with a median age of $\sim 195$ Myr; and (3) region G, a
population of stars with a median age of $\sim$81 Myr. Although Figure 2 shows
that the upper age limit of our sample is 1 Gyr, 90% of our stars fall within
region A. This is evidence that most of the stars we study are around 20 Myr
old. In addition, younger stars are spatially more clustered than older stars,
and stars older than $\sim$100 Myr exhibit a rather smooth distributions
(e.g., Sun et al., 2017b). Therefore, although our UMS sample covers an age
range from a few Myr to just under 1 Gyr, the majority of stars which our
detection algorithm selects is younger than 100 Myr. In Section 3.3 we discuss
age differences across the LMC and how those could affect our analysis.

####  3.1.4 Lifetimes

Figure 5: Structure map of the LMC coloured by the corresponding crossing
times (in Myr). The colour bar is scaled from 5 Myr to 50 Myr, but the values
range from 1 Myr to 160 Myr.

Figure 5 maps each structure labelled by its derived crossing time. Larger
structures have longer crossing times. Allison et al. (2009) suggested that
all stellar structures will be eliminated into the general background density
of a galaxy within a crossing time. Studying the spatial distribution of LMC
stars, Bastian et al. (2009) supported this suggestion. Therefore, there might
be a link between the crossing time of a structure and its lifetime before the
stellar structure dissolves into the LMC field. For comparison, the typical
crossing time in open clusters is close to 1 Myr (Lada & Lada, 2003). The
derived crossing times range from 1 Myr to 160 Myr for the identified young
stellar structures. Thus, they will evolve into a smooth distribution within
$\sim$100 Myr, consistent with what Bastian et al. (2009); Sun et al. (2017b)
derived.

###  3.2 Parameter correlations and distributions

Figure 6: Structure demographics on a logarithmic scale. The diagonals show
the histograms of each parameter. Below the diagonals, bivariate scatter plots
are shown with the respective Pearson correlation coefficients overplotted in
black. The scatter plot data points are labelled by their corresponding
significance level; see the legend on the right.

In this section we show how the parameters we calculated in Section 3.1 are
related. Figure 6 shows pair-wise correlations on a logarithmic scale of the
structures’ radii, numbers of stars, surface mass densities and crossing
times. The diagonal panels show the histograms of each parameter. The lower
plots are bivariate scatter plots with the respective Pearson correlation
coefficients ($\rho$) overplotted. $R$ and $t_{\rm cr}$ have $\rho=1$, from
the definition of the Larson relation. $R$ is strongly correlated with
$N_{*}$, with $\rho=0.94$. $N_{*}$ and $t_{\rm cr}$ have $\rho=0.94$, as
$t_{\rm cr}$ is derived from $R$. The other parameters are not strongly
correlated.

####  3.2.1 Perimeter–Area Relation

Figure 7: Perimeter–area relation of the young stellar structures. The red
solid line is a power-law fit to the structures with area $>1.3\times 10^{3}$
pc2. The vertical black dashed line is drawn at the area threshold of
$1.3\times 10^{3}$ pc-2, which corresponds to a structure with a radius of 20
pc. The blue dashed line shows the perimeter–area relationship that is
consistent with circular structures.

Figure 7 shows the area ($A$) versus perimeter ($P$) relationship of the young
stellar structures. Fitting a power-law slope to the perimeter–area relation
assesses whether the shapes of the structures are more regular (circular) or
more irregular. In a plane, circles follow the relationship $P\propto
A^{0.5}$. For an area $\leq 1.3\times 10^{3}$ pc2, Figure 7 follows the
perimeter–area relationship for circles. We overplot a line with a power-law
slope of $\alpha=0.5$ for points with area $<1.3\times 10^{3}$ pc2 to
demonstrate that they follow this relationship and also that it cuts off the
points above this threshold. The threshold is related to the resolution of the
KDE map (10 pc). We choose $A=1.3\times 10^{3}$ pc2, corresponding to $R=20$
pc, to be greater than the resolution of our KDE map, a locus where Figure 7
is not strongly cut off by the perimeter–area relationship pertaining to
circles. Beyond this area, we fit a power-law slope to the data using least-
squares fitting; we find a best-fitting slope $\alpha=0.72\pm 0.01$. Lovejoy
(1982) followed Mandelbrot’s theory of fractals (Mandelbrot, 1983) to study
the geometry of atmospheric clouds on Earth. They used the formula $P\propto
A^{D_{p}/2}$ to parametrise the perimeter–area of the clouds. In this
equation, $D_{p}=1$ for a circle, increasing towards 2 for more irregular
shapes. Falgarone et al. (1991) applied this relationship to interstellar
molecular clouds. We apply this formula to the slope we fit to Figure 7. The
resulting perimeter–area dimension for the structures is $D_{p}=1.44\pm 0.2$.
This indicates that the structures larger than $R$ = 20 pc have very
irregular, non-circular shapes.

####  3.2.2 Mass–Size Relation

Figure 8: Mass–size relation of the structures. (Top) Scatter plot with
different colours labelled by their significance level (shown in the legend).
(Bottom) Power-law slopes fitted to the points at each significance level,
labelled by their significance levels.

Figure 8 shows the mass ($N_{*}$) versus size ($R$) relation of the young
stellar structures at each significance level. The top panel of Figure 8 is a
scatter plot coloured by each significance level, while the bottom panel shows
power-law slopes fitted to the points through least-squares fitting at each
significance level. The power-law slopes vary from shallow, around
$\kappa=1.49$, at the lowest significance level to steep, around 1.9, at the
highest significance levels. Although the the slope deviates at both ends, it
is very similar from 3$\sigma$ to 8$\sigma$. Fitting a line to the entire
sample, we find a power-law slope, $\kappa=1.50\pm 0.10$.

The change in power-law slope at different significance levels can be
explained using the 2D fractal dimension, $D_{2}$. For masses which are
hierarchically clustered, the mass–size relation inside a fractal follows the
relationship $M\propto R^{D_{2}}$ (Mandelbrot, 1983). The 2D fractal dimension
is related to, but not the same as, the perimeter–area dimension, $D_{p}$ (see
Section 3.2.4). Applying this to a stellar mass–size relation, we can derive
the projected fractal dimension, $D_{2}=\kappa$. Feitzinger & Galinski (1987)
was the first study to employ Mandelbrot’s theory of fractals to stellar
populations by determining the fractal dimensions in young star-forming
regions in galaxies. They also applied their algorithm to artificial
distributions so as to compare their results with those of star-forming
regions. They found the uniform artificial distributions to be rhombic and
quadratic and have fractal dimensions $\sim$2, whereas star-forming regions
have fractal dimensions ranging from $\sim$1.4 to $\sim$1.9. Therefore,
stellar structures with a fractal dimension $\sim$2 have a uniform stellar
surface density, whereas those with smaller fractal dimensions are
characterised by patchy, nonuniform distributions. Applying this framework to
the derived slope of our entire sample, we find $D_{2}=\kappa=1.50\pm 0.10$.

####  3.2.3 Size, Mass and Surface Density Distributions

Figure 9: Size distribution of the structures in bins of 0.1 dex. The
histogram peaks at 12.5 pc and the mean of the entire sample is $\sim$21.4 pc.
The red line is a power-law fit to the data between 13 pc and 700 pc. Error
bars represent Poissonian uncertainties.

Figure 9 shows the size distribution of our young stellar structures in bins
of 0.1 dex. The distribution peaks at $\sim$13 pc. There is a decline in
smaller structures, caused partially by the incompleteness of structures owing
to our detection limit of 10 pc and the criterion that each structure must
have $N_{*}>5$ (see also Sun et al., 2017a, b, 2018). Between 13 pc and 700
pc, the size distribution is well-described by a single power law. We hence
fit a power law to this size range. We will comment in more detail on this
apparent break in the $R$ distribution in Section 4.2.

We fitted the power law using least-squares fitting and found
$\alpha(R)=-1.61\pm 0.20$. We found that changing the bin size of the size
distribution does not significantly affect the derived power-law slope. Using
the power-law slope of Figure 9, we can again derive a 2D fractal dimension.
For a hierarchically structured sample, substructures follow a cumulative size
distribution:

$N(>R)\propto R^{-D_{2}},$ (4)

where $D_{2}$ is the fractal dimension (Mandelbrot, 1983; Sun et al., 2017b,
2018). Equation 4 for a cumulative size distribution is equivalent to a
differential size distribution of

$\frac{{\rm d}N}{{\rm d}\log R}\propto R^{-D_{2}}.$ (5)

Using Equation 5, we derive a fractal dimension $D_{2}=-\alpha(R)=1.61\pm
0.20$. This slope is consistent with that derived form the mass-size relation,
within uncertainties.

As indicated by Figure 9, there are three structures with $R>$ 700 pc.
Considering their boundaries in the context of the LMC, we find that they are
similar in size to the overall galactic structure of the LMC. These larger
structures are most significantly influenced by galactic processes and
therefore do not follow the same power-law size distribution.

Figure 10: $N_{*}$ distribution of the structures in bins of 0.15 dex. The red
line is a power-law fit to the data covering the range from 30 to 30,000
stars. The peak of the histogram is at 21.45 stars. Error bars represent
Poissonian uncertainties.

Figure 10 shows the mass ($N_{*}$) distribution of the young stellar
structures. The mass distribution is also incomplete at smaller values. It
appears from Figure 8 that once $N_{*}>30$ the detection limit of $R=10$ pc
becomes inconsequential. Therefore, we fit a power law from 30 to 30,000
stars. We found that there is no significant change in the fitted power-law
slope of the mass distribution when changing the bin size. When plotting the
logarithmic mass distribution of masses arranged hierarchically in a fractal-
like manner, they should have a power-law slope of $-1$ (Elmegreen &
Falgarone, 1996; Sun et al., 2018). That value is within our margin of error.

There are $\sim$5 structures that have $N_{*}>$ 30,000 and which are not well-
represented by the same power law. This is for the same reason as for the $R$
distribution in Figure 9: these structures are more likely associated with
galactic-scale processes rather than sub-galactic processes.

Figure 11: Surface density distribution of the structures in bins of 0.09 dex
with a lognormal distribution overplotted in red. The mean and standard
deviation of the distribution are included in the legend. Error bars represent
Poissonian uncertainties.

Figure 11 shows the distribution of the surface densities ($\Sigma$) for a bin
size of 0.09 dex. We fit a lognormal distribution to the data, with the mean
and standard deviation listed in the legend. We found that changing the bin
size does not affect the distribution. Simulations of supersonic non-
gravitating turbulent gas in isothermal environments have a lognormal shape in
their column density PDFs (Federrath et al., 2010; Konstandin et al., 2012;
Gouliermis et al., 2017). Therefore, the lognormal shape of Figure 11
indicates that the distribution might be the consequence of turbulence. We
elaborate on this further in Section 4.1.

####  3.2.4 On the Fractal Dimension

Fractal dimensions are used in the literature to show self-similar, irregular
structures of groups of stars and gas. There are different methods of
obtaining a fractal dimension. $D_{p}$ and $D_{2}$ are both types of fractal
dimensions, but they are not equivalent.

As Sun et al. (2017a) noted, $D_{2}$ is the fractal dimension of our stellar
groups’ projections onto the 2D plane perpendicular to the line of sight. A
suggested relation between the 2D and 3D fractal dimensions is $D_{3}=D_{2}+1$
(Beech, 1992). This relation is used in the literature for gas structures, but
it has been shown as improbable for star-forming regions (Beech, 1992; Sánchez
et al., 2005, 2007; Sánchez & Alfaro, 2008). Nonetheless, we find that it is
useful to derive $D_{3}$ for purposes of comparison and to differentiate
between $D_{2}$ and $D_{p}$. Deriving fractal dimensions for gas clouds,
Sánchez et al. (2005); Sánchez & Alfaro (2008) showed that $D_{3}>D_{p}+1$.
Our results agree with their findings, since $D_{p}$ is smaller than both our
derived $D_{2}$ values. Using our two derivations for $D_{2}$ as well as our
derivation of $D_{p}$, we estimate that $2.34\leq D_{3}<2.91$. As Sun et al.
(2017a) noted, typical values in the ISM are $D_{3}=2.4$. Therefore, the
limits of our $D_{3}$ agree with values found in the ISM more generally.

###  3.3 Comparison of young stellar structures in different LMC environments

Figure 12: CMDs of stars in three different LMC environments. All plots
contain the same isochrones, labelled with their age in the leftmost plot.
(Left) CMD of Structure 2, with stars plotted at the 1, 7 and 15 $\sigma$
significance levels. (Middle) CMD of the Bar region, with stars plotted at the
3, 9 and 15 $\sigma$ significance levels. (Right) CMD of the 30 Doradus
region, with stars plotted at the 3, 9 and 15 $\sigma$ significance levels.

In this Section, we inspect three different environments in the LMC: (1) the
Bar region; (2) the 30 Doradus region; and (3) Structure 2, also referred to
as $S_{2}$, which contains the Northwestern Arm and Shapely Constellation III;
(1) and (2) are located in $S_{1}$. $S_{1}$ and $S_{2}$ are the largest and
second largest structures identified, labelled on the right-hand side of
Figure 3.

In Section 3.1.3, we commented on the intrinsic difficulties related to
deriving ages for our young stellar structures. To give an idea of the age
ranges, we plot CMDs of all structures at three different significance levels
in these three regions: see Figure 12. Structure 2 and the 30 Doradus region
both show a locus of stars which peak around the isochrones of 10 Myr and 32
Myr, and which show scatter across the entire selection box. However, the Bar
region shows the peak shifted redwards towards 100 Myr. There is also shown a
cluster of stars around the blue-loop region, corresponding to an age of 100
Myr, in the Bar region. Therefore, the stars in the Bar seem slightly older.
Figure 12 displays an expected result based on star-formation studies of the
LMC. The Bar region has a larger line-of-sight depth, it is more affected by
crowding and contains a tangled mixture of stellar populations. Current and
recent star formation are not concentrated in the LMC Bar (Harris & Zaritsky,
2009; Mazzi et al., 2021). In the LMC, two of the most active star-forming
regions are 30 Doradus and the Constellation III region; for example, star-
formation activity occurring $<$12 Myr ago is dominated by these two regions
(Harris & Zaritsky, 2009). There is some recent star-formation activity in the
Bar, but the most active star-forming periods for the LMC Bar occurred instead
5 Gyr, 500 Myr and 100 Myr ago (Harris & Zaritsky, 2009). This can explain the
concentration of stars around 100 Myr in the middle panel of Figure 12.

Figure 13: Cumulative size distribution of the Bar, 30 Doradus and Structure 2
($S_{2}$) regions in the LMC. A power-law slope is fitted between 15 pc and
200 pc. The power-law slope is the same if we fit all points of the three
regions or if we take the average of the power-law slopes of the three
regions. The error represents the combined fitting error of the three regions.
Error bars are not plotted to avoid confusion.

In order to assess if our results change based on environmental location in
the LMC, particularly since there are significant age differences (relative to
young stars) across the LMC, we ran our analysis separately for each region.
We include the cumulative size distribution of the three regions in Figure 13.
The power-law slope is fitted between 15 pc and 200 pc, and represents the
average of three separate power-law slopes (which agree within 1$\sigma$ with
each other). The error represents the combined fitting error. We also find
that all the power-law slopes and fractal dimensions agree within 1$\sigma$.
Therefore, we show that our results do not change based on LMC environment,
and we are confident that the difference in stellar populations illustrated in
Figure 12 does not affect our conclusions.

##  4 Discussion

###  4.1 Comparisons with studies of the ISM

Gas structures have size scales from 0.1 pc to 1 kpc and exist in galaxies as
a hierarchy of clouds, clumps, cores and filaments (Sánchez & Alfaro, 2008).
Molecular clouds are the birth sites of stars. They show a distinct fractal
nature (Elmegreen & Falgarone, 1996; Blitz & Williams, 1999; Krause et al.,
2020) with larger clouds hosting smaller, more compact clouds. The largest gas
structures are on the same scale as the parent galaxy (Elmegreen et al.,
2000). They have the same spatial scales as our young stellar structures (a
few parsecs to $\geq$1 kpc; see Section 3.1.3).

Gas clouds – including molecular clouds, high-velocity clouds and Hi
distributions – have a perimeter–area dimension from 1.2 to 1.5 (Sánchez et
al., 2007). Our derived perimeter–area dimension is in this range (see Section
3.2.1). Sánchez et al. (2005) found a 3D fractal dimension, $D_{3}$, of
$2.5\leq D_{3}\leq 2.7$ for the Orion A molecular cloud. Considering the Hii
distributions in galactic discs, Sánchez & Alfaro (2008) found 2D fractal
dimensions from 1.5 to 2, corresponding to a 3D fractal dimension of 2.73.
Sánchez et al. (2007) found a 3D fractal dimension of 2.7 in molecular clouds.
These are very similar to the 2D fractal dimensions we derived (see Sections
3.2.2 and 3.2.3) and the 3D fractal dimension ranges we estimated for the
young stellar structures in the LMC (see Section 3.2.4). This is evidence the
LMC young stellar structures inherited their organisation from the fractal
ISM.

On galactic scales, Hi is found to have a fractal nature in the LMC (Kim et
al., 2003) and the SMC (Stanimirovic et al., 1999). Elmegreen et al. (2001)
studied Hi emission in the LMC and suggested that the hierarchical structure
results from supersonic turbulence, based on comparisons with computer
simulations (Elmegreen, 1999; Pichardo et al., 2000). Turbulence is
acknowledged as an important driver of ISM structural evolution. It is
believed to play a major role in the creation of molecular clouds and star
formation (Pingel et al., 2018). Turbulence is created on many spatial scales
and cascades to the smallest scales, as shown by the fractal ISM (Elmegreen et
al., 2001). On larger scales, accretion of the circumgalactic medium and
gravitational instabilities probably trigger turbulence (Krumholz & Burkhart,
2016), whereas on smaller scales stellar feedback in the form of outflows and
supernova explosions probably injects energy (Padoan et al., 2016). In the
SMC, the turbulent properties of Hi are observed to be homogeneous (Nestingen-
Palm et al., 2017), which seems to be owing to large-scale gravitational
driving of turbulence. In contrast, the LMC shows a large diversity in
turbulent properties in Hi (Szotkowski et al., 2019). This reveals regions of
the LMC that are dominated by gravitational driving of turbulence but also
regions on spatial scales smaller than 200–300 pc, where stellar feedback
seems to be strongly influencing the turbulent properties.

The fractal nature of these gas structures has been attributed to the physical
processes that create them in the ISM. Observations of interstellar gas and
dust show a fractal nature in their distributions, which can be compared with
theory to indicate that turbulence is an important process in the ISM.
Simulations of compressively driven supersonic turbulence imply a 3D fractal
dimension of $D_{3}=2.7$ (Federrath et al., 2009). This is in the range of our
derived $D_{3}$ values (see Section 3.2.4) and indicates that the clustering
of star formation in the LMC on scales from 10 pc to $\sim$700 pc is driven by
supersonic turbulence.

It is important to note that hierarchical structure arises from both top-down
(fragmentation of molecular clouds) and bottom-up (injection of energy by
stellar winds and supernovae) processes. Therefore, the smallest structures we
have observed could follow more closely the fractal structure of self-
gravitating molecular clouds, whereas the larger structures could possibly
follow more closely the larger-scale ISM structure dominated by the top end of
the supersonic turbulence spectrum.

Studies of molecular clouds show that those that follow a lognormal surface
density distribution are dominated by supersonic turbulence (Padoan &
Nordlund, 2002; Elmegreen, 2011; Gouliermis et al., 2017). In contrast,
molecular clouds that have a surface density distribution with a power-law
tail are more strongly controlled by self-gravity. We showed that the surface
density distribution of the young stellar structures follows a lognormal
distribution (Section 3.2 and Figure 11). However, there is the caveat that
the absence of a power-law tail in our size distribution could be attributed
to the constraints of our study, namely the resolution of our KDE map of 10
pc. We might not be able to detect the effects of self-gravity in our study,
because we cannot detect the smallest, densest structures. The lognormal shape
of our surface density distribution is a further demonstration that the
congregation of star formation in the LMC is regulated by supersonic
turbulence.

###  4.2 Comparisons with past work on the Magellanic Clouds

Gouliermis et al. (2003) expanded on the stellar detection method introduced
by Gouliermis et al. (2000) and applied it to a $6.5^{\circ}\times
6.5^{\circ}$ area centred on the LMC Bar. With a detection limit of 20 pc,
they detected 494 systems. They separated the systems into bound, intermediate
and unbound using their calculated central densities. Their unbound systems
are described as representing stellar associations. They found a mean diameter
of these of 85 pc. The mean size or peak in the size distribution of a young
stellar system is important, because it could represent a universal scale for
star formation. We will discuss the presence or absence of this characteristic
scale in Section 4.3.

Bastian et al. (2009) studied the spatial distributions of different stellar
populations in the LMC using the optical Magellanic Cloud Photometric Survey
(Zaritsky et al., 1997). They used a fractured minimum spanning tree approach
to break their sample of stars into separate groups. For this, a breaking
distance was applied in which all edges of the tree longer than the breaking
distance were removed. They applied 20 different breaking distances. Their
detection limit was 1.5 pc. From this approach they created a catalogue of
objects compiled from every breaking distance applied and removed duplicates.
Using Two Point Correlation Functions (TPCF) and the $Q$ parameter proposed by
Cartwright & Whitworth (2004), which can differentiate between centrally
concentrated and fractal distributions, they derived a 2D fractal dimension
(see Table 3). From their results, they hypothesised that the stellar
structure evolution in the LMC may be most affected by general galactic
dynamics: the galactic velocity dispersion may erase substructures and render
the stellar distribution uniform. They conjectured that all stellar structures
would be eliminated in a crossing time. From our results, we suggest that the
young stellar structures are not most affected by general galactic dynamics.
It is not the motion across the LMC that erases the young structures. Instead,
it is their internal disequilibrium that causes them to disperse and then
local dynamical effects that smooth the stellar distribution. Additionally,
the LMC’s global velocity dispersion is dominated by older stars which
acquired deviant velocities as a result of dynamical relaxation. Our young
stellar structures do not experience global stirring; they form part of an
orderly, rotating disc (Gaia Collaboration et al., 2021).

Plotting the cumulative size distribution (their figure 1) using four
different breaking distances, Bastian et al. (2009) found no evidence of a
characteristic scale of star formation. They found that the breaking distance
used corresponds to the mean of the distribution in a one-to-one manner.
Therefore, they concluded that there is no evidence of a characteristic scale
of star formation in the LMC and that, in general, characteristic size scales
found in different studies may be the result of the selection of a single
breaking distance or size scale used.

For the cumulative size distribution for their total sample, they note the
existence of an excess at radii $>400$ pc. This break is similar to that seen
in the autocorrelation function of young stars in the LMC (1 kpc; Odekon,
2008) and the scale where a break in the Hi distribution is observed: 100–290
pc (Elmegreen & Elmegreen, 2001; Padoan et al., 2001; Kim & Park, 2007). It
has been suggested that these observed breaks are owing to the disc’s
geometry, for example the scale height of the disc, which is observed to have
a varying thickness from $\sim$50 pc to $\sim$ 400 pc (Szotkowski et al.,
2019). We find a similar ‘break’ in our size distribution (Figure 9) around
300–700 pc. We suggest that this is related to the scale height of the LMC,
the size of the galaxy, and the disc structure. Owing to the break of
$\sim$700 in our size distribution, we infer that around this scale a
transition occurs: below this scale supersonic turbulence dominates, whereas
above this scale global galactic structure dominates.

Sun et al. (2017a, b) studied young stellar structures in two VMC tiles which
contain 30 Dor and the Bar star-forming complex. They found structures ranging
from a few parsecs to 100 pc in size. Sun et al. (2018) studied young stellar
structures in an area of $\sim 28$ square degrees centred on the SMC using 16
VMC tiles. They identified 556 structures on 15 significance levels, whereas
we identified 2846 structures on 15 significance levels. Their structures
range in size from a few parsecs to 1000 pc. We compare our results with Sun
et al. (2017a, b, 2018) in Table 3. We find similar results for the LMC and
SMC. In their size distribution, Sun et al. (2018) found that the peak was at
$\sim$10 pc, similar to the peak we found at $\sim$12 pc. One important
difference is the fitting range they used for their size distribution: from 10
pc to 100 pc. Their upper bound is lower owing to the smaller size of the SMC.
A comparison of our results shows that young star-forming regions with ages
from a few Myr to 100 Myr are similar in the LMC and the SMC in terms of their
properties (size, mass, area, etc.) despite differences in their physical and
chemical environments.

Zivkov et al. (2018) studied 31 PMS structures in the LMC using VMC data.
Fitting a power law to the number distribution of their structures, they
derived a slope of $-0.86\pm 0.12$. This agrees with our value within the
errors.

Table 3 compares the derived parameters $D_{2}$, $D_{p}$, the power-law fit to
$N_{*}$ and the shape of the surface density distribution from this work and
four other similar studies of the Magellanic Clouds. In general, the results
are very similar. We suggest that the differences in $D_{2}$ may be owing to
the method of determination (power-law fit versus the TPCF/$Q$ parameter).
When fitting a power law, the results depend on whether the cumulative or
differential size distribution is considered. Moreover, the slope depends on
the size of the binning scheme and the mathematical fitting algorithm. The
derived exponents for $D_{2}$ are all $<2$, which is expected for a fractal
distribution of stars.

Table 3: Fractal dimensions and power-law slopes for young stars in the Magellanic Clouds. Columns include data about the region of interest and reference and derived $D_{2}$, $D_{p}$, a power-law fit to the number of stars ($\alpha(N_{*})$) and the shape of the surface density distribution ($\Sigma$ Dist.). $D_{2}$ values are further labelled by a superscript which indicates the derivation method: 1mass–size relation; 2size distribution; 3TPCF/$Q$ parameter. Galaxy/Region | $D_{2}$ | $D_{p}$ | $\alpha(N_{*})$ |  $\Sigma$ Dist.
---|---|---|---|---
LMC: 30 Doradus (Sun et al., 2017a) | $1.6\pm 0.3^{2}$ | $--$ | $--$ | $--$
LMC: Bar (Sun et al., 2017b) |  $1.5\pm 0.1^{1}$; $1.5\pm 0.1^{2}$; $1.75\pm 0.01^{3}$ | $--$ | $-1\pm 0.1$ | lognormal
LMC: This work |  $1.50\pm 0.10^{1}$;$1.61\pm 0.20^{2}$ | $1.44\pm 0.20$ | $-0.97\pm 0.18$ | lognormal
SMC: (Sun et al., 2018) |  $1.48\pm 0.03^{1}$;$1.4\pm 0.1^{2}$ | $1.44\pm 0.20$ | $-1\pm 0.1$ | lognormal
LMC: (Bastian et al., 2009) | $\sim 1.8^{3}$ | $--$ | $--$ | $--$

###  4.3 Comparisons with past work on other galaxies

Table 4: Stellar structures detected in other galaxies. Listed are the galaxy and, in parentheses, whether or not the study encompasses the entire galaxy, the distance, the number of stellar structures detected, the resolution of the study (if known), the $R$ distribution peak (converted to radius if indicated originally in units of diameter), $D_{2}$ and (if relevant) the fitting range, $D_{p}$ and the shape of the surface density ($\Sigma$) distribution. The superscripts associated with the $D_{2}$ values indicate the derivation method: 1mass–size relation; 2size distribution; 3autocorrelation function. The subscripts adjacent to galaxy names indicate the reference: 1this work; 2Sun et al. (2018); 3Gouliermis et al. (2015); 4Gouliermis et al. (2017); 5Rodríguez et al. (2016); 6Rodríguez et al. (2018); 7Rodríguez et al. (2019); 8Rodríguez et al. (2020).  Galaxy(-Wide?) | Dist. | Structures | Res. (pc) |  $R$ Peak (pc) | $D_{2}$ | $D_{p}$ |  $\Sigma$ Shape
---|---|---|---|---|---|---|---
LMC1 (yes) | 50 kpc | 2846 | 10 | 13 |  $1.50^{1}$;$1.61^{2}$ (13–700 pc) | $1.44$ | lognormal
SMC2 (yes) | 62 kpc | 556 | 10 | 10 |  $1.48^{1}$;$1.40^{2}$ (10–100 pc) | $1.44$ | lognormal
NGC 65033 (yes) | 5 Mpc | 244 | 80 | 65 |  $1.6^{1}$;$1.5^{2}$ (65–1000 pc);$1.7^{3}$ | $--$ | $--$
NGC 15664 (yes) | 10 Mpc | 890 | 67 | 63 |  $1.54^{1}$;$1.8^{2}$ (63–1000 pc) | $--$ | lognormal, bimodal
NGC 3005 (no) | 1.93 Mpc | 1147 | $--$ | 25 | $--$ | $1.58$ | $--$
NGC 2536 (no) | 3.56 Mpc | 875 | $--$ | 40 | $--$ | $1.50$ | $--$
NGC 2477 (no) | 3.60 Mpc | 339 | $--$ | 55 |  $\sim 1.8^{2}$ (30–150 pc) | $1.58$ | $--$
NGC 24038 (no) | 3.18 Mpc | 573 | $--$ | 40 | $--$ | $1.64$ | $--$

The parameters we will pay most attention to in this section are the fractal
dimensions ($D_{2}$ and $D_{p}$) and the mean/peak (which are closely related)
of the size distributions. As stated in Section 4.2, the mean/peak of the size
distribution in studies of young stellar structures or stellar associations
can be indicative of a universal, characteristic scale of star formation.

First, as regards the characteristic scale (mean/peak of the size
distribution), Mel’Nik & Efremov (1995) studied OB associations within 3 kpc
of the Sun and found a maximum in the size distribution at 15 pc and an
average radius of 20 pc. Bresolin et al. (1998) studied the OB associations of
seven spiral galaxies and found that the size distribution peaked at radii of
20–40 pc. Pietrzyński et al. (2005) compared OB association parameters of 22
galaxies and found a peak of the size distribution at 25–55 pc with an average
of 20–90 pc. Efremov et al. (1987) and Ivanov (1996) found that young stellar
associations have average radii of $\sim$40 pc. Owing to the abundance of
stellar association studies, a radius of $\sim$40 pc has been discussed as a
characteristic galactic scale for star formation (see Gouliermis, 2011).
However, it is important to note that these studies only considered unbound
stellar groups as stellar associations, not smaller, more compact systems in
the star-formation hierarchy.

Gouliermis et al. (2015) and Gouliermis et al. (2017) expanded on the young
stellar structure method used by Gouliermis et al. (2003). The parameters of
their studies are listed in Table 4. They found peaks in their size
distributions of 63–65 pc. This is larger than those found in this paper and
by Sun et al. (2018). We suggest that this is because of the distances to the
galaxies considered and the resolution limit of their studies (80 and 67 pc,
respectively). Their derived $D_{2}$ values are very similar to ours and are
based on a similar fitting range. Gouliermis et al. (2017) found a bimodal
surface density distribution. They noted that the first node is composed
entirely of structures found at the first significance level. This is
different from our unimodal surface density distribution (see Figure 11), but
it did not lead them to reach different conclusions from ours.

Gouliermis et al. (2017) also used the viral theorem to estimate the velocity
dispersions and crossing times of their structures. Because many of their
structures are unbound, their crossing times are upper limits to the true
values, whereas the derived velocity dispersions are lower limits to the true
values. From their 1$\sigma$ to 12$\sigma$ levels, they found an average range
of 0.5–0.8 km s-1 for the velocity dispersion and an average range of 71–406
Myr for the crossing times. The average values we estimated for the velocity
dispersion (see Table 2) are larger by an order of magnitude and display a
different trend: the lowest significance-level structures (on the 1 and
2$\sigma$ levels) have the highest velocity dispersions. We assert that using
the Larson relation leads to more accurate velocity dispersion estimates for
young stellar structures. One reason is that studies show that velocity
dispersions for younger clusters increase as they age (Ramírez-Tannus et al.,
2021). Another is that our values are similar to the velocity dispersion for
RSG stars in the young LMC cluster NGC 2100 (Patrick et al., 2016).

Rodríguez et al. (2016, 2018, 2019, 2020) studied four galaxies in the Local
Group. The parameters of their studies are listed in Table 4. It is important
to note that these four studies are not galaxy-wide studies; they instead
looked at one or more images covering part of the galaxy. They used a path-
linkage criterion technique to detect young stellar groups and create a size
distribution. They also used KDE maps. They found similar $D_{2}$ and $D_{p}$
values to ours and the other studies included in Table 4.

Grasha et al. (2017) studied six local (3–15 Mpc) star-forming galaxies. With
a completeness limit of 1 pc, they studied 3685 young stellar groups. From
application of the TPCF, they found an average $D_{2}$ of 1.17, but they
emphasised that their results vary for each galaxy and that the different
clustering strengths imply different fractal dimensions and different star-
formation hierarchies for each galaxy. They found a characteristic radius from
50 pc to 150 pc. Their average $D_{2}$ is much smaller than those listed in
Tables 3 and 4. This could be owing to the different method they used to
estimate $D_{2}$. The characteristic size is larger than the mean and peak of
our size distribution, but similar to many of those listed in Table 4.

Mondal et al. (2019) studied the dwarf irregular galaxy IC 2574 (3.79 Mpc) and
identified 419 bright far-UV groups. With a resolution limit of 15 pc, they
found that their groups have radii of 15–285 pc and they identified some of
them with stellar associations. Mondal et al. (2021) completed a similar study
for the flocculent spiral galaxy NGC 7793 and identified 2046 young star-
forming clumps. With a detection limit of 6.8 pc, they found radii between 12
pc and 70 pc. These two studies have a similar lower bound to their young
star-forming regions, but a lower upper bound than ours. This is owing to the
sizes of the galaxies considered and the fact that they use far-UV data; in
addition, their definition of a star-forming group is different from our young
stellar structures.

Gusev (2014) studied NGC 628 (7.2 Mpc) and, using a detection method similar
to Gouliermis et al. (2010), detected 297 young stellar structures with a
detection limit of 50 pc. They found a bimodal size distribution, with the
highest peak at a radius of $\sim$30 pc and a lower peak at 300 pc.
Considering the cumulative size distribution, they fitted a power law. They
found a power-law slope of $-1.5$ in three radial ranges: (1) 23–43 pc; (2)
95–300 pc; and (3) 325–450 pc. They stated that they found three
characteristic sizes of young star groups: (1) OB associations with mean radii
of 33 pc; (2) stellar aggregates with mean radii of 120 pc; and (3) stellar
complexes with mean radii of 292 pc. Therefore, they did not find a continuous
power-law slope in the range of 50–1000 pc. In contrast, we found a continuous
power law covering scales of 13–700 pc and a single peak in our size
distribution. These differences could be owing to the detection method, the
different types of distribution used to estimate the power-law slopes or the
fitting algorithm used. No comparable study found similar breaks in their
cumulative or differential size distributions.

Comparing results from the LMC, SMC and the other galaxies referred to in this
subsection, we can see if the hierarchical properties depend on metallicity
and/or environment (e.g., starburst, galaxy interactions, etc.). It seems that
there is no obvious dependence. Therefore, regardless of which physical
processes take place on large galactic scales or which microphysics controls
the small scales, the clustering of star formation between the scale of
clusters and galaxies seems to be universally controlled by turbulence-driven
hierarchical star formation.

##  5 Conclusion

We present the highest-resolution study to date of the observed hierarchical
architecture of groups of young stars in the LMC. This is the first such study
done covering an area of $\sim$120 square degrees across the LMC. The main
results of the study are as follows:

  1. 1.

We identified 2846 young stellar structures on 15 significance levels from a
KDE map (with a resolution of 10 pc) from young stars selected from 68 VMC
tiles.

  2. 2.

We showed the hierarchical organisation of the young stellar structures by
displaying a dendrogram of the second-largest structure.

  3. 3.

We calculated parameters of the young stellar structures: radius, area,
perimeter, surface density, number of stars, crossing time and velocity
dispersion.

  4. 4.

The structures range in size from a few parsecs to 1000 pc, with a mean radius
of $\sim$21 pc.

  5. 5.

The velocity dispersions range from 1.7 to 23 km s-1.

  6. 6.

The crossing times range from 2 to 160 Myr.

  7. 7.

From the perimeter–area relation, we derived a perimeter–area dimension of
$D_{p}=1.44\pm 0.20$.

  8. 8.

From the mass–size relation, we derived a 2D fractal dimension of
$D_{2}=1.50\pm 0.10$.

  9. 9.

The size distribution of the structures has a peak at $\sim$13 pc. From the
size distribution, we derived a 2D fractal dimension of $D_{2}=1.61\pm 20$ in
the range of 13–700 pc.

  10. 10.

From the $N_{*}$ we fitted a power law from 30 to 30,000 stars, with a slope
$\alpha(N_{*})=-0.97\pm 0.18$.

  11. 11.

We showed that the surface density distribution fits that of a lognormal
distribution.

Our distributions and fractal dimensions ($D_{2}$ and $D_{p}$) are similar to
those found previously in the LMC and other galaxies. There are strong
similarities in the fractal dimensions and the surface density distribution to
those of structures in the ISM and numerical simulations. Therefore, we
conclude that our identified stellar structures inherited their fractal nature
from the ISM and that supersonic turbulence is a dominant mechanism in
generating this architecture. Our results support the explanation of
turbulence-driven hierarchical star formation across different galactic
scales. Our findings suggest that star formation in the LMC, on scales from 10
pc to $\sim$700 pc, is scale-free and forms a continuum of star-formation
structures.

## Acknowledgements

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement no. 682115). This research was also supported in part by the
Australian Research Council Centre of Excellence for All Sky Astrophysics in 3
Dimensions (ASTRO 3D), through project number CE170100013. We would like to
thank the Cambridge Astronomy Survey Unit (CASU) and the Wide Field Astronomy
Unit (WFAU) in Edinburgh for providing the necessary data products under the
support of the Science and Technology Facilities Council (STFC) in the U.K.
This study is based on observations made with VISTA at the ESO/La Silla
Paranal Observatory under programme ID 179.B-2003. Finally, this project has
made extensive use of the Python packages: Astropy (The Astropy Collaboration
et al., 2018), matplotlib (Hunter, 2007), NumPy (Harris et al., 2020), pandas
(pandas development team, 2020) and SciPy (Virtanen et al., 2020).

## Data Availability

The entire young stellar structure catalogue is available at
https://github.com/amyelizmiller/VMC_LMC_hierarchical_star_formation. The VMC
images are publicly available, but the PSF VMC photometry used in this work
are propriety to the VMC team, with a data release planned in early 2022; see
https://www.eso.org/sci/publications/announcements/sciann17232.html. The SMASH
data used are part of DR2 and are available at
https://datalab.noao.edu/smash/smash.php#SecondDataRelease.

## References

  * Allison et al. (2009) Allison R. J., Goodwin S. P., Parker R. J., de Grijs R., Portegies Zwart S. F., Kouwenhoven M. B. N., 2009, ApJ, 700, L99
  * Bastian et al. (2007) Bastian N., Ercolano B., Gieles M., Rosolowsky E., Scheepmaker R. A., Gutermuth R., Efremov Y., 2007, MNRAS, 379, 1302
  * Bastian et al. (2009) Bastian N., Gieles M., Ercolano B., Gutermuth R., 2009, MNRAS, 392, 868
  * Bate et al. (1998) Bate M. R., Clarke C. J., McCaughrean M. J., 1998, MNRAS, 297, 1163
  * Beech (1992) Beech M., 1992, Ap&SS, 192, 103
  * Bekki (2017) Bekki K., 2017, MNRAS, 469, 2933
  * Bell et al. (2019) Bell C. P. M., et al., 2019, MNRAS, 489, 3200
  * Bica et al. (2008) Bica E., Bonatto C., Dutra C. M., Santos J. F. C., 2008, MNRAS, 389, 678
  * Blitz & Williams (1999) Blitz L., Williams J. P., 1999, arXiv e-prints, pp astro–ph/9903382
  * Bonatto & Bica (2010) Bonatto C., Bica E., 2010, MNRAS, 403, 996
  * Bresolin et al. (1998) Bresolin F., et al., 1998, AJ, 116, 119
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Bressert et al. (2010) Bressert E., et al., 2010, MNRAS, 409, L54
  * Cartwright & Whitworth (2004) Cartwright A., Whitworth A. P., 2004, MNRAS, 348, 589
  * Cen (2021) Cen R., 2021, ApJ, 906, L4
  * Cioni et al. (2011) Cioni M.-R., et al., 2011, A&A, 527, A116
  * Cross et al. (2012) Cross N. J. G., et al., 2012, A&A, 548, A119
  * Dalton et al. (2006) Dalton G. B., et al., 2006, ] 10.1117/12.670018, pp 62690X–62690X–10
  * Efremov & Elmegreen (1998) Efremov Y. N., Elmegreen B. G., 1998, MNRAS, 299, 588
  * Efremov et al. (1987) Efremov I. N., Ivanov G. R., Nikolov N. S., 1987, Ap&SS, 135, 119
  * El Youssoufi et al. (2021) El Youssoufi D., et al., 2021, MNRAS, 505, 2020
  * Elmegreen (1999) Elmegreen B. G., 1999, ApJ, 527, 266
  * Elmegreen (2000) Elmegreen B. G., 2000, ApJ, 530, 277
  * Elmegreen (2010) Elmegreen B. G., 2010, in de Grijs R., Lépine J. R. D., eds, Vol. 266, Star Clusters: Basic Galactic Building Blocks Throughout Time and Space. pp 3–13 (arXiv:0910.4638), doi:10.1017/S1743921309990809
  * Elmegreen (2011) Elmegreen B. G., 2011, ApJ, 731, 61
  * Elmegreen & Elmegreen (2001) Elmegreen B. G., Elmegreen D. M., 2001, AJ, 121, 1507
  * Elmegreen & Falgarone (1996) Elmegreen B. G., Falgarone E., 1996, ApJ, 471, 816
  * Elmegreen et al. (2000) Elmegreen B. G., Efremov Y., Pudritz R. E., Zinnecker H., 2000, in Mannings V., Boss A. P., Russell S. S., eds, Protostars and Planets IV. p. 179 (arXiv:astro-ph/9903136)
  * Elmegreen et al. (2001) Elmegreen B. G., Kim S., Staveley-Smith L., 2001, ApJ, 548, 749
  * El Youssoufi et al. (2019) El Youssoufi D., et al., 2019, MNRAS, 490, 1076
  * Falgarone et al. (1991) Falgarone E., Phillips T. G., Walker C. K., 1991, ApJ, 378, 186
  * Federrath et al. (2009) Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010, A&A, 512, A81
  * Feitzinger & Galinski (1987) Feitzinger J. V., Galinski T., 1987, A&A, 179, 249
  * Flaugher et al. (2015) Flaugher B., et al., 2015, AJ, 150, 150
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A7
  * Girardi et al. (2005) Girardi L., Groenewegen M. A. T., Hatziminaoglou E., da Costa L., 2005, A&A, 436, 895
  * Gomez et al. (1993) Gomez M., Hartmann L., Kenyon S. J., Hewett R., 1993, AJ, 105, 1927
  * González-Fernández et al. (2018) González-Fernández C., et al., 2018, MNRAS, 474, 5459
  * Gouliermis (2011) Gouliermis D. A., 2011, Phys. Scr, 84, 048401
  * Gouliermis et al. (2000) Gouliermis D., Kontizas M., Korakitis R., Morgan D. H., Kontizas E., Dapergolas A., 2000, AJ, 119, 1737
  * Gouliermis et al. (2003) Gouliermis D., Kontizas M., Kontizas E., Korakitis R., 2003, A&A, 405, 111
  * Gouliermis et al. (2010) Gouliermis D. A., Schmeja S., Klessen R. S., de Blok W. J. G., Walter F., 2010, ApJ, 725, 1717
  * Gouliermis et al. (2015) Gouliermis D. A., et al., 2015, MNRAS, 452, 3508
  * Gouliermis et al. (2017) Gouliermis D. A., et al., 2017, MNRAS, 468, 509
  * Grasha et al. (2017) Grasha K., et al., 2017, ApJ, 840, 113
  * Grudić et al. (2018) Grudić M. Y., Guszejnov D., Hopkins P. F., Lamberts A., Boylan-Kolchin M., Murray N., Schmitz D., 2018, MNRAS, 481, 688
  * Gusev (2014) Gusev A. S., 2014, MNRAS, 442, 3711
  * Ha et al. (2021) Ha T., Li Y., Xu S., Kounkel M., Li H., 2021, ApJ, 907, L40
  * Harris & Zaritsky (2009) Harris J., Zaritsky D., 2009, The Astronomical Journal, 138, 1243
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hodge (1973) Hodge P. W., 1973, AJ, 78, 807
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Irwin et al. (2004) Irwin M. J., et al., 2004, ] 10.1117/12.551449, pp 411–422
  * Ivanov (1996) Ivanov G. R., 1996, A&A, 305, 708
  * Keller et al. (2001) Keller S. C., Grebel E. K., Miller G. J., Yoss K. M., 2001, The Astronomical Journal, 122, 248–256
  * Kim & Park (2007) Kim S., Park C., 2007, ApJ, 663, 244
  * Kim et al. (2003) Kim S., Staveley-Smith L., Dopita M. A., Sault R. J., Freeman K. C., Lee Y., Chu Y.-H., 2003, ApJS, 148, 473
  * Konstandin et al. (2012) Konstandin L., Girichidis P., Federrath C., Klessen R. S., 2012, ApJ, 761, 149
  * Krause et al. (2020) Krause M. G. H., et al., 2020, arXiv
  * Krumholz & Burkhart (2016) Krumholz M. R., Burkhart B., 2016, MNRAS, 458, 1671
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Larson (1995) Larson R. B., 1995, MNRAS, 272, 213
  * Lovejoy (1982) Lovejoy S., 1982, Science, 216, 185
  * Lucke & Hodge (1970) Lucke P. B., Hodge P. W., 1970, AJ, 75, 171
  * Mandelbrot (1983) Mandelbrot B. B., 1983, The fractal geometry of nature / Benoit B. Mandelbrot, updated and augmented [ed.] edn. W.H. Freeman New York
  * Mazzi et al. (2021) Mazzi A., et al., 2021, MNRAS, 508, 245
  * Mel’Nik & Efremov (1995) Mel’Nik A. M., Efremov Y. N., 1995, Astronomy Letters, 21, 10
  * Mondal et al. (2019) Mondal C., Subramaniam A., George K., 2019, AJ, 158, 229
  * Mondal et al. (2021) Mondal C., Subramaniam A., George K., Postma J. E., Subramanian S., Barway S., 2021, ApJ, 909, 203
  * Nestingen-Palm et al. (2017) Nestingen-Palm D., Stanimirović S., González-Casanova D. F., Babler B., Jameson K., Bolatto A., 2017, ApJ, 845, 53
  * Nidever et al. (2017) Nidever D. L., et al., 2017, AJ, 154, 199
  * Nidever et al. (2021) Nidever D. L., et al., 2021, AJ, 161, 74
  * Odekon (2008) Odekon M. C., 2008, ApJ, 681, 1248
  * Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
  * Padoan et al. (2001) Padoan P., Kim S., Goodman A., Staveley-Smith L., 2001, ApJ, 555, L33
  * Padoan et al. (2016) Padoan P., Pan L., Haugbølle T., Nordlund Å., 2016, ApJ, 822, 11
  * Patrick et al. (2016) Patrick L. R., Evans C. J., Davies B., Kudritzki R. P., Hénault-Brunet V., Bastian N., Lapenna E., Bergemann M., 2016, MNRAS, 458, 3968
  * Pichardo et al. (2000) Pichardo B., Vázquez-Semadeni E., Gazol A., Passot T., Ballesteros-Paredes J., 2000, ApJ, 532, 353
  * Pietrzyński et al. (2005) Pietrzyński G., Ulaczyk K., Gieren W., Bresolin F., Kudritzki R. P., 2005, A&A, 440, 783
  * Pingel et al. (2018) Pingel N. M., Lee M.-Y., Burkhart B., Stanimirović S., 2018, ApJ, 856, 136
  * Ramírez-Tannus et al. (2021) Ramírez-Tannus M. C., et al., 2021, A&A, 645, L10
  * Rodríguez et al. (2016) Rodríguez M. J., Baume G., Feinstein C., 2016, A&A, 594, A34
  * Rodríguez et al. (2018) Rodríguez M. J., Baume G., Feinstein C., 2018, MNRAS, 479, 961
  * Rodríguez et al. (2019) Rodríguez M. J., Baume G., Feinstein C., 2019, A&A, 626, A35
  * Rodríguez et al. (2020) Rodríguez M. J., Baume G., Feinstein C., 2020, A&A, 644, A101
  * Rolleston et al. (2002) Rolleston W. R. J., Trundle C., Dufton P. L., 2002, A&A, 396, 53
  * Rubele et al. (2012) Rubele S., et al., 2012, A&A, 537, A106
  * Rubele et al. (2015) Rubele S., et al., 2015, MNRAS, 449, 639
  * Sánchez & Alfaro (2008) Sánchez N., Alfaro E. J., 2008, ApJS, 178, 1
  * Sánchez et al. (2005) Sánchez N., Alfaro E. J., Pérez E., 2005, ApJ, 625, 849
  * Sánchez et al. (2007) Sánchez N., Alfaro E. J., Pérez E., 2007, ApJ, 656, 222
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Skowron et al. (2021) Skowron D. M., et al., 2021, ApJS, 252, 23
  * Stanimirovic et al. (1999) Stanimirovic S., Staveley-Smith L., Dickey J. M., Sault R. J., Snowden S. L., 1999, MNRAS, 302, 417
  * Subramanian & Subramaniam (2009) Subramanian S., Subramaniam A., 2009, A&A, 496, 399
  * Sun et al. (2017a) Sun N.-C., et al., 2017a, ApJ, 835, 171
  * Sun et al. (2017b) Sun N.-C., et al., 2017b, ApJ, 849, 149
  * Sun et al. (2018) Sun N.-C., et al., 2018, ApJ, 858, 31
  * Sun et al. (2020) Sun N.-C., Maund J. R., Crowther P. A., 2020, MNRAS, 497, 5118
  * Sutherland et al. (2015) Sutherland W., et al., 2015, A&A, 575, A25
  * Szotkowski et al. (2019) Szotkowski S., et al., 2019, ApJ, 887, 111
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration Price-Whelan A. M., Sipőcz B. M., Günther H. M., et al., 2018, AJ, 156, 123
  * Urrutia Zapata et al. (2019) Urrutia Zapata F., Fellhauer M., Alarcón Jara A. G., Matus Carrillo D. R., Aravena C. A., 2019, MNRAS, 489, 2746
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vázquez-Semadeni et al. (2016) Vázquez-Semadeni E., González-Samaniego A., Colín P., 2016, MNRAS, p. stw3229
  * Wan et al. (2020) Wan Z., Guglielmo M., Lewis G. F., Mackey D., Ibata R. A., 2020, MNRAS, 492, 782
  * Zaritsky et al. (1997) Zaritsky D., Harris J., Thompson I., 1997, AJ, 114, 1002
  * Zivkov et al. (2018) Zivkov V., et al., 2018, A&A, 620, A143
  * de Grijs et al. (2014) de Grijs R., Wicker J. E., Bono G., 2014, AJ, 147, 122
  * pandas development team (2020) pandas development team T., 2020, pandas-dev/pandas: Pandas, doi:10.5281/zenodo.3509134, https://doi.org/10.5281/zenodo.3509134

##  Appendix A Filling in the Gap with SMASH data

Figure 14: SMASH CMD in the rectangular region around the gap of VMC data. The
black contours show the density of the cross-matches of UMS stars. The red
polygon is the region we used to select stars. Figure 15: Pixel-to-pixel
comparison of the pixels of the SMASH and VMC KDE maps. The points represent
mean $\mu_{\text{VMC}}$ values in 15 equally spaced $\mu_{\text{SMASH}}$ bins.

Figure 1 (in Section 2.1) shows an absence of data centred at ($\Delta$RA,
$\Delta$Dec) = $(0.65,-2.48)$ with a width of 0.108∘ and a height of 1.445∘.
To fill this gap, we followed the method of Sun et al. (2018). They filled a
gap in VMC data of the SMC using $UBV$ photometric data. We filled the gap in
the LMC coverage using optical data from SMASH Data Release 2 Nidever et al.
(2021). SMASH data were gathered with the Dark Energy Camera (DECam; Flaugher
et al., 2015) and cover the Magellanic System in fields of $\sim$3 deg2 in the
$ugriz$ passbands. A complete description of the survey strategy and data
reduction can be found in Nidever et al. (2017, 2021). To use SMASH data, we
first defined a rectangular region around the gap with dimensions of
2${}^{\circ}\times 1.26^{\circ}$. This region encloses a region of 0.1 pc to
0.2 pc around the gap; in this region there are 1119 UMS stars selected from
the VMC data around the gap. We identified this region of sky in the SMASH
survey. To select point-like sources from the SMASH data, we made cuts using
three morphology parameters: $\tt{CHI}<3$, $|\tt{SHARP}|<1$ and
$\tt{PROB}>0.8$. Using this selection, we cross-matched with UMS stars in the
same rectangular region. There were 1094 matches in this region. Figure 14
shows a CMD of SMASH data in our rectangular region, with the black contours
representing the density of cross-matches and the red polygon showing our
selection of UMS stars. There are 1937 SMASH stars in this region of CMD
space.

Figure 16: KDE maps of the rectangular region centred at ($\Delta$RA,
$\Delta$Dec) = $(0.65,-2.48)$ with size 2${}^{\circ}\times 1.26^{\circ}$.
(Left) KDE map of the VMC data. (Middle) KDE map of the calibrated SMASH data.
(Right) Combined map of the VMC and SMASH data.

We constructed surface density maps of the rectangular region using both VMC
and SMASH data. The left-hand panel of Figure 16 shows the VMC KDE map of this
region. We calibrated the pixels of the SMASH KDE map, $\mu_{\text{SMASH}}$,
by doing a pixel-to-pixel comparison, shown in Figure 15. We fitted a line to
this comparison and found $\mu_{\text{VMC}}=0.65\mu_{\text{SMASH}}$. The
middle panel of Figure 16 shows the calibrated SMASH KDE map. To combine these
two maps, we used equation (2) of Sun et al. (2018):
$\mu_{\text{comb}}=\mu_{\text{VMC}}\times w+\mu_{\text{SMASH}}\times(1-w)$,
where $w$ is a weight assigned to each point depending on position: $w=0$ if
the pixel falls within 10 pc from the gap edge, increasing by 0.01 every 10 pc
from the gap edge until $w=1$. The result is shown in the right-hand panel of
Figure 16. Finally, we merged the combined KDE map with our full VMC surface
density map, shown in the top panel of Figure 3. The resulting map has a mean
stellar density of $3.21\times 10^{-3}$ pc-2, a median stellar density of
$7.58\times 10^{-5}$ pc-2 and a standard deviation of $7.65\times 10^{-3}$
pc-2.

##  Appendix B Comparison with other surveys

Figure 17: (Left) Our final structure map in greyscale with three surveys of
LMC young stellar systems overplotted. (Right) Our structure map before the
final selection criterion is applied, with the same three surveys overplotted.

We compare our structures with published galaxy-wide catalogues of young
stellar systems in the LMC. If the provided coordinates fall within the
1$\sigma$ isodensity contours, we estimate that we detected the system. If
not, we estimate that system is not detected in our study.

The first survey used for this comparison is the Lucke & Hodge (1970)
catalogue, containing 122 associations. They defined stellar associations as
loose groups of bright blue stars, which they identified from $B$ and $V$
photographic plates. We detect 95% of their associations in our study.

The next survey was published by Gouliermis et al. (2003). They identified
young stellar systems from a photographic $U$ plate centred on the LMC. Their
detection method is similar to ours: they identified contours of young systems
on significance levels if they are 3$\sigma$ above the local mean stellar
density. Then they classified the systems based on their calculated central
densities. They provide two catalogues, one of their unbound and one of their
intermediate stellar systems, which they describe as associations and open
clusters, respectively. To compare with our structures, we impose that their
stellar systems have at least five stars. After applying this criterion, their
open cluster catalogue contains 231 systems and their association catalogue
has 148 systems. We detect 77% of the open cluster catalogue and 65% of the
association catalogue.

The final catalogue we compare with is the Bica et al. (2008) catalogue, which
is the most complete catalogue of stellar systems in the Magellanic Clouds.
Their catalogue contains extended objects which were identified by eye on
ESO/SERC $R$ and $J$ Sky Survey Schmidt films, as well as catalogues of
systems available in the literature at the time of publication. We compare
with objects classified as ‘CN’ or ‘cluster in nebula’, ‘NC’ or ‘nebular with
probable embedded cluster’, and ‘NA’ or ‘nebular with embedded association’.
We chose these three classes because they were classified as the youngest
stellar systems in the LMC by Bonatto & Bica (2010), with reported ages of
$\leq$10 Myr. There are 81 CN systems and we detect 95%. There are 183 NC
systems and we detect 93%. There are 817 NA systems and we detect 88%.

From the above analysis, it transpires that we have the best overlap with the
Bica et al. (2008) and Lucke & Hodge (1970) catalogues. The Lucke & Hodge
(1970) and Gouliermis et al. (2003) criteria are most similar to those used in
our study, so we provide our structure map with the boundaries of their
systems overplotted in Figure 17. On the left we plot our final structure map,
and on the right we plot our structure map before the final selection
criterion is applied. The final selection criterion is that all of our
structures on the 1$\sigma$ and 2$\sigma$ levels must contain at least one
structure at the 3$\sigma$ and higher levels. This means that all of our
structures on the 1$\sigma$ and 2$\sigma$ levels have at least 10 stars.
Although Gouliermis et al. (2003) used a very similar detection algorithm, the
slight algorithm differences and differences in data caused them to find more
structures in certain low-density regions. As Figure 17 shows, our algorithm
detects most structures in the Gouliermis et al. (2003) catalogue, but after
we apply our final selection criteria we exclude $\geq$30 of them. We find
that it is necessary to apply the final selection criterion, because we
identify structures at levels above the median stellar density, which is
effectively zero. In contrast, Gouliermis et al. (2003) identified structures
above the local mean. We must make a compromise to cull spurious objects and
we made the decision to apply global selection criteria. Furthermore, the
purpose of our study is to analyse the fractal properties, not to identify new
stellar clusters or associations.

Finally, we estimated how many structures are newly identified objects by
cross-matching our catalogue with the catalogues mentioned in this section.
The central coordinates of our structures correspond to the geometric
centroid. If our structure does not match within 0.5 arcsec the coordinates in
one of those other surveys, we estimate that the structure of interest is a
newly discovered stellar system. Doing this, we estimate that 45% of our
structures are newly discovered systems.