© The Authors 2025

1 Introduction

The formation of globular clusters (GCs) remains an open question, with no definitive explanation agreed upon within the community at present. They are dense and compact stellar systems with no residual gas, characterized by typical stellar surface densities of 10⁴ M_⊙ pc⁻² and typical sizes of a few parsecs (Baumgardt et al. 2020; Baumgardt & Vasiliev 2021). A considerable fraction of GCs presents several chemical anomalies and/or multi-modal metallicity distributions (Gratton et al. 2004; Carretta et al. 2009) that are nowadays interpreted as a signature of multiple populations (Gratton et al. 2012, Bastian & Lardo 2018, Milone & Marino 2022 and references therein). Also, the general consensus is that they lack dark matter (e.g., Moore 1996; Conroy et al. 2011). To form stellar clusters with these high densities and compact sizes, it is reasonable to speculate that specific conditions must be met. For instance, it is necessary to accumulate and convert large amounts of gas into stars efficiently within a short timescale and a confined volume, which requires environments of exceptionally high pressure and density (pressure P/k > 10⁷ K cm⁻³ and density n > 10⁴−10⁵ cm⁻³; Elmegreen & Efremov 1997; Kruijssen 2014, 2015; Renaud et al. 2015). However, despite understanding these necessary conditions, many other aspects of GC formation remain unresolved. Key uncertainties include the specific mechanisms driving the formation of multiple stellar populations (Bastian & Lardo 2018) and the question of whether GCs initially formed within dark matter halos (Peebles 1984; Boley et al. 2009; Trenti et al. 2015; Kimm et al. 2016; Renaud 2018) that were later stripped by tidal interactions (Mackey & van den Bergh 2005) or if they formed in environments devoid of dark matter (Naoz & Narayan 2014; Phipps et al. 2020).

The majority of observed GCs is relatively old, with ages exceeding 10 Gyr (Krauss & Chaboyer 2003; Forbes & Bridges 2010), indicating that they must have formed in a young Universe, prior to reionization, where gas-rich and dense conditions were prevalent (Tacconi et al. 2013; Wisnioski et al. 2015; Dessauges-Zavadsky et al. 2019). In this context, high-redshift observations are therefore essential, providing critical insights into the formation and evolution processes behind these objects. Indeed, data from the Hubble Space Telescope (HST) and, more recently, the James Webb Space Telescope (JWST) have enabled the discovery of numerous compact and massive sources in lensed fields at very high redshifts (Vanzella et al. 2017, 2019; Calura et al. 2021; Meštrić et al. 2022; Mowla et al. 2022; Vanzella et al. 2023; Claeyssens et al. 2023; Adamo et al. 2024; Messa et al. 2024). These sources exhibit extreme properties, with stellar densities up to 10⁵ M_⊙ pc⁻², or even 10⁶ M_⊙ pc⁻², confined to small physical sizes ranging from 2–3 pc to a few tens of parsecs. While the exact nature of these dense systems remains uncertain, their similarities to present-day GCs suggest the possibility of an evolutionary link, implying that they could serve as progenitors of local GCs.

To enhance our understanding of these compact sources and their connection to GCs, as well as to complement the wealth of observational data available, it is crucial to develop a comprehensive theoretical framework that accurately captures the physical processes responsible for their formation but also considers the broader cosmological context, essential to tracing their evolution. In this respect, cosmological simulations are vital to bridge this gap in knowledge. However, to effectively characterize the small-scale dynamics occurring on parsec scales, they must achieve exceptionally high spatial and mass resolutions (Kimm et al. 2016; Ma et al. 2020) and must incorporate advanced feedback models (Agertz et al. 2013; Dale 2015; Ceverino et al. 2017; Weinberger et al. 2017; Pillepich et al. 2018; Hopkins et al. 2018; Marinacci et al. 2019).

Historically, in the context of modeling small structures such as dwarf galaxies, most simulations have represented star clusters as single particles, operating at spatial resolutions of approximately 10 pc and mass resolutions of 10³ M_⊙. Such values may fail to adequately sample these small systems, leading to an inadequate modeling of the subparsec processes involved in star formation and the resulting feedback. To accurately capture their complexity, simulations must achieve sub-pc spatial resolutions, employ stellar particle masses of at most 10 M_⊙, or (ideally) model individual stars that directly sample the initial mass function (IMF; e.g., Sormani et al. 2017) to get a proper stochastic representation of star formation (Andersson et al. 2020; Calura et al. 2022; Deng et al. 2024). However, fully resolving star clusters in early galaxies remains a challenge, primarily due to the dramatically increased computational costs associated with achieving such high resolutions across entire cosmic volumes.

Recent studies have advanced our understanding of the formation of GCs and their connection to high-redshift stellar clumps. As an example, by means of cosmological simulations, Ma et al. (2020), Sameie et al. (2023), and Garcia et al. (2023) have suggested that bound star clusters arise in high-pressure and high-density environments, highlighting the importance of the star formation model and how it influences cluster formation. Similarly, Calura et al. (2022) and Calura et al. (2025, hereafter referred to as C24) conducted zoom-in cosmological simulations with a subpc resolution, examining the effects of various stellar feedback models on the formation of high-redshift clumps. The connection of GCs and dark matter has been addressed by Gutcke (2024, as well as Kimm et al. 2016) demonstrating that multiple star formation bursts occur in low-mass dark matter halos and that dark matter can be efficiently lost by stripping, making these systems resemble GC systems. However, despite these improvements, many limitations are still present. Most simulations are halted at high redshift (z > 5–6) due to the high computational cost (e.g., Kimm et al. 2016), fail to replicate the observed densities and/or sizes, cover incorrect mass ranges, or a combination of thereof (e.g., Ma et al. 2020). The underlying reasons for these discrepancies are not entirely understood and may be the result of various factors, such as resolution limitations, missing physical processes, or processes that have been inadequately implemented. Additionally, even when dealing with simulations that can successfully generate dense stellar aggregates with compact sizes and correct masses, they frequently struggle to accurately reproduce other dynamical or chemical signatures (Kimm et al. 2016).

In this work, we focus on the detailed analysis of a specific simulation within the framework of the SImulating the Environment where Globular clusters Emerged (SIEGE) project (Calura et al. 2022). This project encompasses a series of cosmological, hydrodynamical simulations at a subparsec resolution, specifically designed to investigate the physical conditions that lead to the formation of dense stellar clusters. The high resolution of these simulations is required to model feedback processes at the scale of individual stars, which is crucial for accurately reproducing the complex dynamics that govern early star cluster formation.

The simulation analyzed in this work follows the evolution of a small cosmological volume centered around a massive dark matter halo, which is the progenitor of a dwarf galaxy, down to a redshift of z = 6.14. At this redshift, the halo has a virial mass of approximately 4 × 10⁹ M_⊙ (or, equivalently, 4 × 10¹⁰ M_⊙ within three virial radii). The primary aim of this study is to investigate the structural properties of the stellar agglomerates formed within the volume and to assess how representative these clusters are of the dense stellar clumps observed in high-redshift lensed fields. To address this, we will analyze the scaling relations and structural properties of the clusters formed in the simulation and track their evolution as a function of redshift. As presented in C24, the feedback model implemented here features a high star formation efficiency, which strongly promotes the formation of compact and dense star clusters. This characteristic directly influences the compactness and density of the resulting systems (C24).

This paper is structured as follows. In Section 2, we provide detailed information about the simulation parameters and setup. Section 3 describes the algorithms employed to identify the stellar clusters, while Section 4 outlines the methods used to extract the structural parameters of these systems. In Section 5, we present our results and in Section 6, we discuss them in the context of the existing literature. In Section 7, we present our conclusions.

2 Simulation setup

The simulation examined in this study was first introduced by C24 within the SIEGE framework. They carried out and presented a series of cosmological simulations aimed at investigating how varying stellar feedback models influence the properties of the earliest star clusters formed at high redshift. They specifically examined the effects of different stellar winds models, varying star formation efficiencies, and distinct IMFs. Below, we summarize the main characteristics of the simulations presented by C24, with a focus on the one that is most relevant to our study. For further details regarding the simulation setup, additional information can be found in Calura et al. (2022), Pascale et al. (2023), and C24.

All simulations from C24 were conducted using the adaptive mesh refinement (AMR) hydrodynamic code RAMSES (Teyssier 2002) and evolved down to z = 10.5, with the exception of the one considered in this work, evolved down to a redshift of z = 6.14. The simulations assume a flat ΛCDM cosmological model with a matter density of Ω_m = 0.276 and a Hubble constant of H₀ = 70.3 km s⁻¹ Mpc⁻¹ (Sharov & Vorontsova 2014; Omori et al. 2019), a convention we will maintain throughout our analysis. The simulation box encompasses a L = 5 Mpc h⁻¹ comoving volume. The dark matter mass resolution is 165 M_⊙ per particle in the zoom-in region, resulting in a total of approximately 2 × 10⁸ dark matter particles. The use of the AMR scheme combined with a relatively small simulation box enables subparsec resolution, with a maximum refinement level of 21¹, corresponding to a spatial resolution of 0.3 h⁻¹ pc at a redshift of 6.14.

As discussed in Calura et al. (2022), all simulations of the SIEGE project were designed to reproduce the properties of the star-forming complex observed in the lensed field strongly magnified by the galaxy cluster MACS J0416.1–2403 (the Cosmic Archipelago; Vanzella et al. 2019; Calura et al. 2021; Messa et al. 2025; Vanzella et al. 2024) at z = 6.14. The complex comprises of several nucleated (<10 pc) stellar agglomerates, with masses from a few 10⁵ M_⊙ up to 5 × 10⁶ M_⊙, with a total mass of no more than a few 10⁷ M_⊙. On the basis of existing stellar-to-halo mass relations for low-mass halos (Ma et al. 2019, see also Fig. 2 of Pascale et al. 2022), a stellar system of this mass is thus expected to be hosted by a dark matter halo with a virial mass of a few times 10⁹ M_⊙, as in the simulation.

Figure 1 shows a representation of the evolution with the redshift of the cosmic volume captured in the simulation. The top, middle, and bottom panels illustrate, the central part of the simulation box at redshifts of 10.5, 8, and 7, respectively. Throughout the panels, the filamentary structure of dark matter is clearly discernible. As the dark matter collapses, it forms increasingly massive knots that are likely to become sites of star formation. Each panel also includes the projected stellar density distribution of the relevant region for a comparison with dark matter. The most massive halos visible at high redshifts begin to form stars at their centers, while being connected by dark matter filaments that gradually increase in size. The two massive halos visible in the top panel, linked by a dark matter bridge, move towards each other and nearly merge by a redshift of 8 (middle panel). By a redshift of 7, these halos have almost completely merged and are about to accrete a third system forming the ≃10⁹ M_⊙ virial mass halo at z = 6.14.

It is worth noting that despite different implementations of the feedback model, the large-scale properties of dark matter remain consistent across the various simulations of C24. Thus, although the dark matter density distributions in the figure are derived from the simulation analyzed in this study, they are representative of the general behavior across all simulations. Significant differences are observed instead in the small-scale properties of the stellar component (in terms of star formation histories, densities, and size of the formed stellar aggregates; C24).

The physical processes incorporated in the simulations include atomic radiative cooling from hydrogen, helium, and metals in photoionization equilibrium with a redshift-dependent ionizing ultraviolet (UV) background (Haardt & Madau 1996), assuming reionisation starting at z = 8.5. The simulations implement the formation of individual stars with the direct IMF sampling method (Sormani et al. 2017). Stars are formed at a rate of $${\dot{\rho }}_{\star }=\frac{{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}}{{\tau }_{\mathrm{S}\mathrm{F}}},$$(1)

with ρ_gas the gas density of a cell and τ_SF the timescale of star formation. The timescale of star formation is parametrized as ${\tau }_{\mathrm{S}\mathrm{F}}=\frac{{\tau }_{\mathrm{F}\mathrm{F}}}{{ϵ}_{\mathrm{S}\mathrm{F}\mathrm{E}}}$, with ${\tau }_{\mathrm{F}\mathrm{F}}\equiv \sqrt{\frac{3\pi }{32G{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}}}$ as the free fall time, ϵ_SFE as the star formation efficiency per free fall time, and G is the gravitational constant. Gas cells eligible for star formation are those with a temperature below 2 × 10⁴ K and that simultaneously reach a density threshold that depends on redshift. The threshold is set based on the requirement that no more than 90% of gas in a cell is converted into stars, with this amount further limited to be, at most, 32 M_⊙. Given a physical resolution of 0.3 pc h⁻¹ at z = 6.14, the condition translates into a density threshold of 10⁴ particles cm⁻³ (Yaghoobi et al. 2022; C24). At such densities and scales, the gas that is typically available for star formation is only enough for a few stars, which motivates the need for individual star formation. Since stars are formed individually, the simulations incorporate feedback from individual stars in the form of stellar winds (SWs) and supernovae (SNe). Also, to prevent artificial radiative losses of the energy injected by the stars, a delayed cooling mechanism tends is empleyed (Teyssier et al. 2013), ensuring an avoidance of numerical overcooling and associated over-abundant star formation.

The simulations in C24 are four and they differ in the implementation of specific aspects of the star formation and/or stellar winds model. Following the nomenclature established in C24, these simulations are (all sharing the same initial conditions):

FC22: The FC22 simulation investigates the impact of stellar winds from massive stars implemented as an impulsive release of energy and mass into the interstellar medium (ISM) at the time of the formation of a stellar particle. The IMF is sampled according to a Kroupa (2001), with a star formation efficiency of ϵ_SFE = 0.1.

WINDS: Compared to the FC22 simulation, the WINDS simulation modifies the stellar winds model by adopting a continuous injection of energy into the surrounding gas throughout the lifetime of the massive stars. The star formation efficiency is ϵ_SFE = 0.1 and the IMF is a Kroupa.

THIMF: In the THIMF simulation, the stellar winds model is consistent with that used in the WINDS simulation; however, the IMF employed is top-heavy, described by $$\varphi (m)=\frac{\mathrm{d}N}{\mathrm{d}m}=\text{\hspace{0.17em}constant.\hspace{0.17em}}$$(2)

This choice is motivated by studies indicating a potential over-abundance of massive stars in early stellar populations (Marks et al. 2012; Charbonnel et al. 2014; Denissenkov & Hartwick 2014; Cameron et al. 2024).

SFE1: In the SFE1 simulation, the IMF follows a Kroupa, and the stellar winds model is consistent with that used in the WINDS simulation. However, this simulation adopts a star formation efficiency of ϵ_SFE = 1, in contrast to the ϵ_SFE = 0.1 assumed in the other simulations (see Section 6.1).

Our focus is on the SFE1 simulation because, according to C24, it consistently generates a significant number of realistic high-density, compact stellar clusters. While C24 offers a comparative analysis of all four simulations, our analysis specifically examines SFE1 due to its apparent relevance in reproducing the characteristics of these stellar agglomerates. This is qualitatively shown in Fig. 2, which presents a selection of star-forming regions from this simulation. Each panel focuses on a different redshift and such clusters are clearly evident as numerous small, high-density knots scattered around the relevant regions. All these knots have sizes of a few parsecs and can reach surface densities of several 10³ M_⊙ pc⁻². A more thorough look at into Fig. 2 is left to Section 3.

Fig. 1 Evolution with redshift of the projected dark matter density within the high-resolution simulation domain. Top, middle, and bottom panels show the projected dark matter density distribution in the central region of the simulation box at a redshift of z = 10.5, 8, and 7, respectively. In all panels, brighter colors indicate higher-density regions. The stellar projected density maps corresponding to the same regions are overlaid for comparison, allowing for a visual examination of the relationship between dark matter and stellar structures across different cosmic epochs. The stellar clumps formed have been highlighted with black circles.

Fig. 2 Examples of stellar density distributions of star-forming regions taken at three different redshifts. The top-right, top-left, and lower maps correspond to z = 10.5,7, and 6.14, respectively. The star-forming systems shown in the top-right and top-left panels correspond to the stellar counterparts of the most massive halos visible in the top and bottom panels of Fig. 1. In each panel, the small insets display the star formation histories of the corresponding regions. All clumps identified within the time windows color-coded in the small insets are marked with circles and displayed with the same colors on the corresponding projected stellar density map.

3 Identification of stellar clumps

A key aspect of this work is the accurate and efficient identification of small stellar agglomerates across the large simulated volume and across redshift. To this end, we developed a method specifically tailored to this simulation. We focused all of the subsequent analysis on eight representative redshifts, namely, z = 15.5, 13, 10.5, 10, 9, 8, 7, and 6.14. For each of the considered redshifts, the identification method was carried out in two distinct steps. First, the simulation volume was divided into non-overlapping regions by identifying the major areas where stars have formed. In the second step, the stellar particles belonging to each of these regions (tracing massive dark matter halos) were grouped into age bins that were independently scanned to search for stellar clumps. Both steps make extensive use of the Hierarchical Density-Based Spatial Clustering of Applications with Noise (Campello et al. 2013) algorithm, as implemented in the software library HDBSCAN (McInnes et al. 2017).

For any given n-dimensional dataset, HDBSCAN constructs a minimum spanning tree (MST) from the distance-weighted graph of the data points, once an appropriate metric has been assumed (Euclidean in our case). An MST is a tree that connects all data points (nodes) with the shortest possible total distance (edge weights) without creating any loops. After constructing the MST, HDBSCAN varies the density threshold to build a hierarchy of clusters. This allows the algorithm to detect clusters of different densities, making it particularly effective for datasets with complex and noisy structures. HDBSCAN is an extension of the Density-Based Spatial Clustering of Applications with Noise (DBSCAN; Ester et al. 1996; Schubert et al. 2017), and one of its advantages is that it does not require us to specify a linking length (ϵ); this is, namely, a characteristic distance used to evaluate whether elements of the dataset belong to the same cluster, but it evaluates stable clusters based on how long they persist as the density threshold changes. Persistent clusters are considered features of the dataset.

Delving into the details of our clump identification process, in the first step, HDBSCAN was run on a randomly selected subsample of stellar particles. Using a subsample is essential to reducing computational costs, especially at lower redshifts, when the number of stellar particles becomes extremely large (>10⁷ at z < 8). Here, we detected clusters directly operating on the MST setting a relatively large density threshold. A large density threshold in the MST allows us to extract only a few, well-populated clusters. Once a region is located, we determined its center as the median particles position and we defined a characteristic radius as the largest Cartesian distance from the region’s center to any of its particles. For cases where the identified regions overlap or are concentric, we handle them as follows. We let r and R be the radii of two overlapping regions, with r < R, and d be the distance between their centers; (i) if d < R − r, the smaller region is eliminated, as it is entirely within the larger one; (ii) if R − r ≤ d < R, the radius of the larger region is expanded to R + r to include the smaller one which was thus removed; (iii) if R ≤ d < R + r, the two regions were merged into one. The new center was calculated as a weighted average of the original centers (with weights proportional to the cubes of the original radii). The new radius is set to the distance between the original centers.

As discussed in detail in the following sections, a key feature of the ϵ_SFE = 1 simulation is that star formation occurs in short, intense bursts. Each burst typically leads to the formation of a distinct stellar agglomerate. Thus, to enhance the effectiveness of our clump-finding algorithm, we leveraged this characteristic by incorporating the ages of star particles into the clustering process. In the second step, for each identified region, we further grouped the stars into age bins of fixed width and ran HDBSCAN separately on the particle positions within each age bin.

HDBSCAN requires us to specify one key parameter, namely: the expected number of objects per cluster. In our implementation, this parameter is set to 500, which enables the identification of smaller-mass agglomerates. Given that stellar particles cannot be less massive than 0.5 M_⊙, this setting allows for the identification of stellar agglomerates with a minimum mass of 250 M_⊙. When applying HDBSCAN to the different age bins, each run typically identifies a large number of clusters. However, not all identified clusters correspond to genuine small, compact stellar agglomerates; some may represent stellar streams, isolated groups of particles, or background noise. Also, the same clump can, in principle, appear in different age bins if it experiences continuous or multimodal star formation. To accurately identify true agglomerates and avoid multiple detections, we employed the following refinement process. We first computed the center of mass for each identified cluster with an iterative scheme (shrinking sphere method; Power et al. 2003). Specifically, we first get an initial estimate of the center of mass using only the stars classified as members by HDBSCAN within the relevant age bin. Subsequently, a refined center of mass is determined by including all particles located within a 5 pc radius of the previous center of mass estimate, regardless of their classification as members or age bin. We then assigned an average density computed using the 100 nearest particles to this mass center, without any age distinction. In our framework, a cluster is thus identified in space by a single point representing its center of mass. Clusters are only retained if their average density exceeds a threshold of 25 M_⊙ pc⁻³. The threshold of 25 M_⊙ pc⁻³ is determined empirically: lower values tend to classify non-genuine agglomerates as clusters, as verified by visual inspection, while higher values tend to exclude progressively denser structures. As a final step, we performed a cross-match of the centers of mass within the star-forming region of interest to identify possible multiple detections of the same clump. If a cross-match was detected, such as when the same stellar structure undergoes two or more star formation bursts, we eliminated one of the two detections.

The clustering algorithm was applied independently to all eight redshifts included in this study. As a result, clumps identified at redshift z_i are not guaranteed to be identified in subsequent redshifts z_i+1 < z_i. Also, we set the age bin width to 25 Myr, but we tested the clustering algorithm with a narrower bin width of 12.5 Myr, obtaining convergent results regarding the number of identified agglomerates. After completing the entire identification procedure, we obtained a set of candidate clumps for each redshift. In the following section, we evaluate the structural properties of these clumps (see Section 4) and conduct a second inspection based on these properties to identify and eliminate any potential remaining false detections.

Dividing the simulation box into star-forming regions and age bins offers several advantages. First, the algorithm runs approximately five times faster compared to executing it directly on all stellar particles in the simulation (≃480 CPU hours vs. ≃80 CPU hours at z = 6.14). Second, it is significantly more efficient in terms of memory usage, allowing the algorithm to operate effectively without requiring high-memory computing clusters. Lastly, segmenting by age bins enhances cluster identification, leading to a higher number of distinct clusters. Without this age-based division, it can be challenging to differentiate among clusters, particularly when they appear as overdensities in areas of low particle density.

Figure 2 shows a selection of stellar projected density distributions of three star-forming regions taken at different redshifts. The stellar map shown in the bottom panel is extracted from z = 6.14 and it depicts the stellar system populating the central and most massive halo at this redshift. The top two panels instead correspond to z = 10.5 (right) and 7 (left), illustrating the same stellar system as in the bottom panel, but during different evolutionary stages, before merging with other massive progenitors. Each star-forming region is populated by a multitude of small, dense knots, each extending over a scale of a few parsecs. Although it is not shown here for the sake of simplicity, these knots are present across the entire simulation, throughout every star-forming region, with maximum densities reaching nearly 10⁴ M_⊙ pc⁻². The abundance of these small and dense agglomerates underscores the necessity for a highly effective clump-finding algorithm. The small panels attached to each projected density map display the star formation history of the corresponding region. The three vertical bands represent age bins, each spanning 25 Myr. All clumps identified within these time windows are highlighted with circles and shown in different colors on the corresponding projected stellar density map. As anticipated, the star formation history is highly fragmented, characterized by frequent bursts of star formation, as seen in the numerous peaks within these panels.

4 Determination of stellar structural parameters

The mass, size, and any derived quantity (e.g., average surface density) of each identified clump can be determined by fitting the clump’s binned spatial density distribution to a well-motivated model. The model surface density profile is parametrized by $${\mathrm{\Sigma }}_{\star }(R)=\mathrm{\Sigma }(R)+{\mathrm{\Sigma }}_{\mathrm{B}\mathrm{G}},$$(3)

where Σ_BG is a constant background added to a Sérsic (1968) surface density model $$\mathrm{\Sigma }(R)={\mathrm{\Sigma }}_0\exp [-b(n){\left(\frac{R}{R_{\mathrm{e}}}\right)}^{1/n}],$$(4)

with $${\mathrm{\Sigma }}_0=\frac{b^{2n}}{2\pi n\mathrm{\Gamma }(2n)}\frac{M_{\star }}{R_{\mathrm{e}}^2}$$(5)

and (Ciotti & Bertin 1999) $$b(n)\simeq 2n-\frac{1}{3}+\frac{4}{405n}+\frac{46}{25515n^2}.$$(6)

In the previous equations, Γ is the Gamma function, M_⋆ is the total stellar mass of the clump, n is the Sersic index (a measure of the system’s concentration), and R_e is the effective radius (i.e., the radius that contains half of the projected mass), such that $$2\pi {\int }_0^{R_{\mathrm{e}}}\mathrm{\Sigma }(R)RdR=\frac{M_{\star }}{2}.$$(7)

For each clump, we considered three independent projections, along the canonical x, y, and z axes. For each projection, the fit is conducted as follows. First, we determine the number of radial bins (N_bins) as the square root of the number of particles within a 10 pc radius from the clump center. The bins span the radial range [0.1, 10] pc, evenly spaced in logarithm. Using this binning scheme, we can compute the clump’s projected density profile, which is fitted to get a first estimate of its effective radius. Using the derived R_e, we adjusted the radial bins to cover the range of [0.1, 10]R_e, we determine a new density profile and subsequently fit it again. The fit was performed by minimizing the figure of merit: $${\chi }^2=\sum _{i=1}^{N_{\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{s}}}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{o}\mathrm{b}\mathrm{s},i}-{\mathrm{\Sigma }}_{\star }\left(R_i\right)}{\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{o}\mathrm{b}\mathrm{s},i}}\right)}^2.$$(8)

Thus, for each projection of each clump, the total number of free parameter is four: n, R_e, M_⋆, and Σ_BG. The triplet {R_i, Σ_obs,i, ΔΣ_obs,i} comprises, respectively, the average distance of the i-th bin, its surface density and uncertainty, and the sum in Equation (8) extends over the total number of bins, with i = 1, ..., N_bins. The minimization of Equation (8) was performed using the Levenberg-Marquardt algorithm (Levenberg 1944; Marquardt 1963).

Following the determination of these structural parameters, we conducted a thorough and additional inspection of the clumps to further refine the sample by eliminating potential outliers. The exclusion criteria include clumps for which the Sérsic fit is unsuccessful in at least one of the three projections, those exhibiting half-mass radii estimates that deviate by more than a factor of 2 in one of the projections compared to the average, and any agglomerates misclassified as a large background region by visual inspection. These additional filtering steps refine the results and eliminates outliers. After applying these criteria, we obtained a sample of clumps, summarized in Table 1.

After this process, each clump is assigned three triplets of total mass, effective radius, and Sérsic index, corresponding to fits from the three different projections. Since these quantities are inherently singular (e.g., a clump should have one mass, not three), we compute an average by taking the mean of the three fitted values, with error bars representing their dispersion². The same method is applied to other derived quantities, such as the average surface density of the clump.

Figure 3 presents a selection of clumps at three distinct redshifts: z = 15.5, 13, and 8, with each row corresponding to one of these redshifts. For each redshift, we chose two clumps, a less massive one shown on the left and a more massive one on the right. Each clump is shown in three different projections, along with the results of the Sérsic fit. Also, in the small insets in each panel we show the corresponding projected stellar density distributions from which the profiles have been computed. Apart from a few cases, the clumps presented in the figure are broadly representative of the entire clump population and do not exhibit significant variations. The majority of clumps have been very well-fitted by a Sérsic profile, indicating a consistent spatial structure (see also Section 5.2). Additionally, as can be appreciated from the projected density maps, the most massive clumps in particular tend to maintain an approximately spherical shape, without significant deviations in their morphology.

To assign the additional characteristic property of age to the identified clumps, we considered all particles within twice the average half-mass radius calculated across different projections and determined the age as the mass-weighted average of the individual particle ages. This method ensures that the computed age accurately reflects the dominant stellar population within the clump. As we demonstrate in detail below, most clumps exhibit minimal age dispersion and have sharply peaked star formation histories (SFHs), rendering a very well defined concept of age.

Fig. 3 Clump projected density distributions. Top panels: projected stellar density profiles (black squares with errobars) for two clumps selected at z = 15.5. Each set of three panels shows the projected density along different lines of sight: along the z, y, and x axes from left to right, respectively (the x z plane is the plane of Fig. 1). The blue curves represent the best-fitting Sérsic + background models, while the vertical lines indicate the position of the effective radius, Re. In the middle and right panels of each triplet we added (with a grey line) the Sérsic + background model shown in the left panel to highlight differences among projections. The curves and the binned profiles shown are those obtained at the end of the two step fitting procedure (see Section 4). The small insets display the projected density map along the same direction as in the main panel, with the blue circle centered at center of the clump and having a radius equal to Re. The clump in the left panel has (M⋆, Re, n) = (1.25 × 104 M⊙, 0.67 pc, 1.24), while the one in the right panel has (M⋆, Re, n) = (1.79 × 104 M⊙, 0.76 pc, 1.47). Middle panels: same as the top panels, but for clumps at z = 10.5. The clumps have (M⋆, Re, n) = (2.88 × 103 M⊙, 1.00 pc, 1.63) – left – and (M⋆, Re, n) = (3.04 × 104, 2.13 pc, 2.69) – right. Bottom panels: Same as the top and middle panels, but for clumps at z = 8, with (M⋆, Re, n) = (3.03 × 103 M⊙, 1.57 pc, 1.62), and (M⋆, Re, n) = (1.83 × 104 M⊙, 1.71 pc, 2.27) on the right. In each panel, the orange shaded area shows the simulation resolution (Δx) at that redshift, Δx ≡ 5h−1/[2l(1 + z)], with l = 21 and h = 0.703.

5 Results

Here, we present the key results of our work. We begin investigating whether the detected clumps are gravitationally bound and whether they form in dark matter halos (Section 5.1). We analyze their structural properties (Section 5.2), study the evolution of various scaling relations with redshift, including the clump mass function (CMF, Section 5.3) and the size–mass relation (Section 5.4), and address the problem of clumps dissolution (Section 5.5).

5.1 Dynamical state of stellar clumps

Although the clustering algorithm is designed to detect dense and compact stellar systems, it does not inherently guarantee that these objects are gravitationally bound, as they may simply represent unbound stellar associations, especially right after their formation. To further investigate their dynamical state, we evaluate the parameter Π, an empirical measure of boundedness introduced by Gieles & Portegies Zwart (2011), defined as $$\mathrm{\Pi }\equiv \frac{t_{\text{age\hspace{0.17em}}}}{t_{\text{cross\hspace{0.17em}}}},$$(9)

where t_cross is the crossing time and t_age is the age of the clump. Π reflects the relationship between a clump’s dynamical time and its age, with values larger than 1 suggesting bound systems, and lower values indicating unbound or transient associations. The crossing time is estimated as $$t_{\text{cross\hspace{0.17em}}}=10\sqrt{\frac{R_{\mathrm{e}}^3}{G{M}_{\star }}}.$$(10)

The middle-left panels of Fig. 4 display the distribution of the Π parameter at three representative redshifts: z = 13, 9 and 7, illustrating the evolution of the clumps’ dynamical state across different epochs. In this case, we limited our analysis only to clumps younger than 50 Myr for two main reasons: (i) in Gieles & Portegies Zwart (2011) the Π parameter is calibrated only for objects below this age threshold; (ii) clumps older than 50 Myr can reasonably be considered bound, as unbound systems would have dissolved on much shorter timescales. The values of Π is almost always larger than 1 throughout the entire redshift range, indicating that the vast majority of these stellar clumps is not transient, unbound structures but rather form as gravitationally bound systems and maintain this state throughout their evolution.

The issue of boundedness becomes more intriguing when considering dark matter. The question of whether GCs and massive high-redshift clumps actually form within dark matter halos remains a subject of debate. While the presence of dark matter may aid in forming bound GC-like objects (Kimm et al. 2016), the mechanism by which they would subsequently lose their dark matter to resemble present day GCs poses some challenges. Current theories suggest that if formation within dark matter halos is the underlying scenario, tidal stripping through interactions with other systems is the primary channel for dark matter loss (e.g., Gutcke 2024).

Given that the Π parameter in Equation (9) measures boundedness primarily for self-gravitating, single-component systems, it is essential to consider how and if dark matter affects our understanding of these clumps. To assess the contribution of dark matter, we estimated the ratio M_dm (3R_e)/M_⋆ for each clump and redshift, where M_dm (3R_e) is the dark matter mass within three times the clumps’ effective radius. The distribution of this ratio at the same three redshifts as for the distribution of Π is shown in the top right panels in Fig. 4. Across all redshifts and clumps, the ratio consistently remains low, rarely exceeding 0.5. The choice of 3R_e as the characteristic scale for this measurement was motivated by the need to account for dark matter distribution while mitigating issues related to the discrete representation of dark matter particles in the simulation. This result underscores the minimal role of dark matter within the clumps, thereby justifying the use of Equation (9) to assess their bound status. The distributions in the top-right panels of Fig. 4 show that the clumps are not dark matter-dominated at the redshifts considered. However, they do not provide insight into whether the clumps initially formed within substantial dark matter halos and later lost dark matter due to processes such as stripping. To investigate this effect we show the relation between the ratio M_dm (3R_e)M_⋆ and the clump ages in the top-left triplet of panels in Fig. 4. If clumps were formed in dark matter-dominated halos and later experienced dark matter loss, we would expect a correlation where younger clumps exhibit higher ratios. However, as shown in the figure, no such correlation is evident at any redshift. This result strongly supports the conclusion that the clumps are currently in environments with minimal dark matter content and formed in regions that were already dark matter-poor.

This result is crucial as it demonstrates the simulation’s ability to form relatively massive stellar systems, reaching masses as high as 10⁵ M_⊙−10⁶ M_⊙ (see Section 5.2), without requiring substantial dark matter halos for formation and stability over time. This aspect will be investigated further in a future work (Pascale et al., in prep.).

5.2 Structural properties of stellar clumps

The remaining triplets of panels in Fig. 4 offer distributions of key structural parameters of stellar clumps across redshifts, based on the Sérsic profile fitting of Section 4. The panels display the distributions of the Sérsic index n (middle right panel), effective radius R_e (bottom left panel), and average surface density (bottom right panel), defined as ${\mathrm{\Sigma }}_{\mathrm{e}}\equiv \frac{M_{\star }}{2\pi R_{\mathrm{e}}^2}$. The middle right panels reveal no significant evolution in the distribution of the Sérsic index n with redshift. The distribution maintains a similar shape, with the index stable and ranging from $n={1.64}_{-0.54}^{+0.64}$ at z = 13 to $n={1.44}_{-0.44}^{+0.59}$ at z = 7, where the uncertainties represent the 16th and 84th percentiles of the distributions. These values span a wide range, indicating profiles with varying degrees of curvature. Small indices suggest cored profiles in the central regions and truncated in the outer parts, while large indices imply more centrally concentrated and shallower outer profiles.

While the distribution of the Sérsic index keeps constant across redshifts, the size distribution – bottom left panels – displays substantial evolution. At z = 13, the distribution is peaked at smaller radii, with a median effective radius of 1 pc, indicating that clumps formed at early times tend to be more compact. As redshift decreases the distribution shifts towards larger radii. The median value of the R_e distribution increases to 3 pc, with some clumps reaching sizes of up to approximately 8 pc.

The shift towards larger effective radii coincides with a decrease in surface density, as shown by the average surface density distributions shown in the bottom right panel of Fig. 4. Indeed, while at high redshifts clumps exhibit very high average surface densities, with a median distribution of 460 M_⊙ pc⁻² and few clumps approaching 6000 M_⊙ pc⁻², as the redshift decreases, there is a notable shift towards lower surface densities. By z = 7, the average surface density drops to 100 M_⊙ pc⁻², with some values falling to as low as 30 M_⊙ pc⁻². This trend, along with the evolution of effective radii, suggests that stellar clumps appear larger over time, regardless of whether they form anew or evolve from pre-existing structures. This expansion may be attributed to two main factors: on one hand, the dynamic evolution, such as the hierarchical merging of smaller structures into larger objects; on the other hand, it could reflect a redshift dependence of gas conditions that promotes the formation of larger and less dense structures at lower redshifts. This aspect will be examined in greater detail in Section 5.5. Finally, we point out that additional processes not included in our simulations such as mass loss due to stellar evolution would cause the stellar systems to undergo an initial expansion.

In the context of the study of the structural properties of stellar clumps, it is important to emphasize that the spatial resolution of our simulation (from ≃0.2 pc at z = 15.5 to ≃0.5 pc at z = 6.14) is close to the characteristic size of the smallest stellar clumps formed. While a formal convergence test would ideally require a higher-resolution run, such simulations are currently computationally prohibitive. We also note that lowering the resolution would not offer a meaningful test of robustness, as the reduced gravitational accuracy would artificially inflate clump sizes, preventing a reliable characterization of their structure. In this regard, the sizes we report should be interpreted as upper limits: resolving smaller spatial scales could potentially reveal even more compact and denser stellar clumps.

Fig. 4 Clump parameters across redshifts. The triplets of panels illustrate: the relation between the dark-matter-to-stellar-mass ratio and clump ages (top-left), the distributions of the dark matter mass within three effective radii relative to the total stellar mass, Mdm (3Re)/M⋆ (top-right), the distributions of the boundedness parameter Π for clumps younger than 50 Myr (middle-left), the distributions of the effective radius Re (middle right), the distributions of the Sérsic index n (bottom left), and the distributions of the average surface density ${\mathrm{\Sigma }}_{\mathrm{e}}\equiv \frac{M_{\star }}{2\pi R_{\mathrm{e}}^2}$ (bottom-right). Each panel within a triplet represents distributions at different redshifts, with z = 13 shown in blue (top), z = 9 in purple (middle), and z = 7 in orange (bottom). The colored bands in the top left triplet of panels highlight clumps younger than 50 Myr. These clumps are used to build the log Π distributions in the middle-left triplet of panels. For clarity, only three of the eight redshifts analyzed in this work are displayed. The vertical dashed line in all triplets (except for the top left) show the median value of the corresponding distribution.

5.3 Clump mass function

A more quantitative analysis of the dependencies between structural parameters of clumps as a function of redshift is presented in Fig. 5. Here, the two sets of panels display, from top to bottom, the CMF, the size–mass relation, and the relationship between size and Sérsic index, while each column is a different redshift. We start focusing on the CMF shown at the top with blue histograms. In each panel, we report with a black dashed line the $\frac{\mathrm{d}N}{\mathrm{d}{M}_{\star }}\propto M_{\star }^{-2}$ power-law relation, while the orange line indicates the best fit to the CMF in the relevant redshift bin, used to quantify changes in the CMF. We do not attempt to fit the CMF at z = 15.5 since only four stellar clumps have formed at this redshift. The slope of the best-fit line is reported in the legend. It is generally expected that the CMF has a power-law behavior, with a typical slope of −2, reflecting the scale-free nature of gas fragmentation of structure formation (Elmegreen & Efremov 1997; Guszejnov et al. 2018). In our case the CMF keeps a slope consistently around −2 within the errorbars, ranging from −2.03 at redshift 13, down to redshift 7, where it is ≃−2.20. We also show in Fig. 5 the CMF of young clumps, selected with an age criterion of t_age < 50 Myr, to make it comparable with observations, where clumps are identified through emission from young stars (e.g., UV emission). At the highest redshift (z = 13), the CMF of young and all clumps almost coincide by definition, as all clumps are inherently young. As redshift decreases, the CMF of young clumps is also stable around a slope of ≃ − 2, in agreement with observations of local (Zhang & Fall 1999; Hunter et al. 2003; de Grijs et al. 2003; Rodruck et al. 2023) or higher redshift young massive clumps (Dessauges-Zavadsky & Adamo 2018; Livermore et al. 2012). The largest differences between the slopes of young and global CMFs are observed at the lowest redshifts, where the global CMF is increasingly dominated by older clumps, particularly at lower masses. Here (z = 7 and 6.14), the CMF of young clumps also deviates more than 1σ from −2, with a difference of 0.3–0.4. The corresponding slope values are provided in Table 1.

Unlike previous studies that observe CMF slope variations and link them to fluctuations in star formation activity (Garcia et al. 2023), we do not observe a significant trend in our simulation. In Fig. 6 we show the SFH of the entire simulation box (top panel) alongside the cumulative mass formed as a function of time (bottom panel). The SFH is computed using bins of varying sizes, with the black curve representing larger bin widths to emphasize the overall mean trend rather than its burstyness, as highlighted by the red line. Up to approximately z = 11, the SFH remains low, averaging a few times 10⁻² M_⊙ yr⁻¹, a value modest compared to the higher star formation rates at lower redshifts. After z = 10.5, the average star formation rate increases by almost a factor of 10, reaching typical values of 0.1 M_⊙ yr⁻¹ and showing a fluctuating pattern with several bumps and sporadic short bursts of ≃0.5 M_⊙ yr⁻¹. For instance, the increase in star formation activity at t = 0.75 Gyr (z = 7.3) corresponds to the major merger of the two massive central stellar systems shown in the middle panel of Fig. 1. This merger event ultimately forms the large system depicted in the top left panel of Fig. 2.

At a redshift of 6.14, we excluded from the fit to the CMF the massive system with a stellar mass of 10⁶ M_⊙ visible in the bottom right panel of Fig. 5. The formation of this massive system is also responsible for the peak in the SFR of 10 M_⊙ yr⁻¹ visible in Fig. 6 at redshift 6.8. This system exhibits singular characteristics, such as an unusually high mass and small size. Also, as highlighted in Section 5.1, this clump forms without a significant contribution from dark matter. Given its peculiar properties and its striking structural and kinematic similarities to local GCs, a dedicated analysis of this system is addressed in a separate study (Pascale et al., in prep.).

5.4 Size–mass relation versus redshifts

The central rows of panels to the top and bottom of Fig. 5 show the clumps size–mass relation as a function of redshift. A correlation between mass and size is evident, with smaller clumps populating the lower-mass region and larger clumps occupying the higher-mass region, although we consistently observe a large scatter in the relation. To quantify this dependence, we have included two curves in each panel: a linear fit in the log R_e − log M_⋆ plane, shown with an orange dashed line, and a line connecting the median R_e computed in each mass bin, represented by a solid red curve. The linear fit follows: $$\log R_{\mathrm{e}}=m_{R_{\mathrm{e}\mathrm{f}\mathrm{f}}-M_{\star }}\log M_{\star }+q_{R_{\mathrm{e}\mathrm{f}\mathrm{f}}-M_{\star }},$$(11)

with m_Reff–M⋆ and q_Reff–M⋆, which are the slope and normalization of the relation, respectively. As shown in Fig. 4, at lower redshifts, the region containing larger sizes becomes more populated, suggesting the formation of increasingly massive structures. Consequently, also the clump size–mass relation exhibits a mild evolution with redshift: after a first increase of the slope in the redshift range z = 15.5–10.5, where it evolves from 0.3 to 0.7, it stabilizes and it keeps approximately the constant value of 0.5 down to redshift 6.14. The values of m_Reff–M⋆ and q_Reff–M⋆ are reported in Table 1 for each redshift bin considered. As already noted in Section 5.2, the spatial resolution of the simulation is close to the characteristic size of the smallest clumps identified. As a result, being able to resolve even smaller spatial scales may potentially lead to a steepening of the relation, increasing the population of low-mass, small-sized clumps.

Although Fig. 4 (top right panel) indicated that the distribution of the Sérsic stays relatively constant across redshifts, the bottom rows of panels in Fig. 5 highlight a weak correlation, with a significant amount of scatter, between the Sérsic index and R_e. In this case we did not provide any fit, but, to highlight this dependence, we show the median Sérsic index per mass bin. Smaller clumps tend to have a median Sérsic index close to 1, which suggests an exponential profile, while larger and more massive clumps display a higher median Sérsic index, with also less scatter.

5.5 Clump dynamics and dissolution

In a simulation where stars are treated as a collisionless component, there are fundamental limitations that prevent the accurate representation of certain physical phenomena. These limitations, by design, encompass collisional processes, which are essential when studying the long-term, dynamical evolution of dense star clusters and GCs. Being able to somewhat account for these phenomena is especially important in simulations like ours, where the dynamics of individual stars are resolved and a high resolution is achieved.

Two-body relaxation is a critical process in the dynamical evolution of star clusters (see e.g., Spitzer 1987, Heggie & Hut 2003). In densely packed environments, stars frequently engage in gravitational interactions, leading to a diffusion process that can cause stars to escape the cluster’s gravitational influence, a phenomenon leading to a decrease in the total mass of a GCs and in some cases to its eventual dissolution. Several factors can accelerate a cluster’s evaporation and dissolution, such as the effects of the tidal truncation due to the external tidal field of the host galaxy, tidal shocks due the time variation of this tidal field, and stellar escape due to the cluster’s expansion associated with mass loss due to stellar evolution and primordial gas expulsion (see e.g., Heggie & Hut 2003, for a review of various dynamical processes contributing to a cluster’s evaporation). A detailed investigation of the evolution of the clumps formed in our simulations is beyond the scope of this paper and will be presented in a future paper; here, we only present a few general considerations and estimates on the dynamics and dissolution of these systems.

The results of a number of studies (see e.g., Chernoff & Weinberg 1990; Vesperini & Heggie 1997; Vesperini 1997; Baumgardt & Makino 2003; Gieles & Baumgardt 2008) have shown that mass loss due to two-body relaxation and the dissolution timescale depend mainly on the cluster’s mass and the strength of the tidal field of the host galaxy. A cluster’s dissolution may be further accelerated by early mass loss due to stellar evolution in particular for low-concentration clusters (see e.g., Chernoff & Weinberg 1990; Fukushige & Heggie 1995).

The complex and strong time variation in the tidal field in the high-redshift external environment (as in the case of one of our simulations) can also significantly accelerate the mass loss and dissolution of clusters. An example of the dynamical effects associated with the high-redshift environment was presented in Li & Gnedin (2019). On the basis of a semi-analytical calculation of the effects of these early time variation of the tidal field, that study showed that most clusters either suffer a significant mass loss during these early evolutionary phases or undergo complete dissolution. In particular their analysis reveals that none of the clusters with masses smaller than 10⁵ M_⊙ survive until the present-day. Most of the clumps formed in our simulation have masses smaller than ~2 × 10⁴ M_⊙ and, although a detailed study of their dynamics will be addressed in a separate paper, we expect the vast majority of them to dissolve before z = 0. Even without considering the effects of the time variation of the high-redshift tidal field, we can derive a very approximate estimate of the dissolution time of these clumps by using the expression of the dissolution timescale obtained by Baumgardt & Makino (2003), as follows: $$\frac{t_{\text{diss\hspace{0.17em}}}}{\mathrm{M}\mathrm{y}\mathrm{r}}\equiv 1.03{\left[\frac{N}{\ln (0.02N)}\right]}^{0.82}\frac{R_G}{\mathrm{k}\mathrm{p}\mathrm{c}}{\left(\frac{V_G}{220\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^{-1}.$$(12)

In the above equation, R_G represents the distance of the clump from the center of its host, V_G is the circular velocity at that location, and N is the total number of stars within the clump. If we calculate, just as an example, the dissolution time by adopting a mild tidal field equivalent to that of the Milky Way at R_G = 6 kpc (and V_G = 220 km s⁻¹), we find that all (or at least the vast majority) clumps dissolve before z = 0: the fraction of disrupted clusters range from 100 percent at z = 15.5 and to 99.5 percent at z = 6.14³. Interestingly, at a redshift of 6.14, the sole clump with a dissolution timescale significantly exceeds 13 Gyr (already highlighted in Section 5.3) and is characterized by a mass of 10⁶ M_⊙ and distinct properties that set it apart from the rest of the population.

We strongly emphasize that this calculation provides just an approximate estimate of the timescale of clump dissolution but a more detailed modeling of the clump dynamics including additional processes, such as the tidal shocks associated with the time variation of the external tidal field, are to further accelerate the clump dissolution and confirm the large fraction of disrupted clusters before z = 0 (e.g., Li & Gnedin 2019).

Finally, we conclude this section with a discussion of the possible role of clump mergers in the evolution of the clump population. Merging of smaller stellar systems into more massive clusters represents another critical process that not only affects the low-mass end but also may contribute to the enhancement over time of the high-mass end of the CMF. If mergers occur on timescales shorter than the dissolution time, they may further reduce the overall impact of dissolution on the CMF, despite its potential influence on low-mass populations of clusters. To quantify this aspect, we turn to Fig. 7, where we present a size–mass relation at z = 7, similar to that shown in Fig. 5, but now color-coded according to the clumps’ ages. Blue indicates younger clumps, while red represents older clumps, with ages computed as described in Section 4. The dashed isocontours reveal a clear trend with age: younger clumps tend to occupy regions of larger sizes and masses, while older clumps are found in areas with smaller sizes and masses. This age dependence indicates that the formation of more massive and larger structures is not primarily driven by merging of smaller clumps, which would populate the top-right part of the relation with older clumps. This is especially true in light of the typical star formation histories of the clumps. The middle panel of Fig. 6 shows the star formation history of a selection of the clumps identified at z = 7, with each color indicating a different clump. The SFHs of these clumps reveals sharp, short-lived episodes of star formation, with each burst lasting only a fraction of a Myr. These SFHs generally display single peaks, rather than multiple bursts. To better highlight this aspect, the lower panels of Fig. 6 present the distribution of log σ_SF for all selected clumps at three representative redshifts. Here, σ_SF is defined as $${\sigma }_{\mathrm{S}\mathrm{F}}^2=\frac{\int \mathrm{S}\mathrm{F}\mathrm{H}(t){(t-t_{\mathrm{a}\mathrm{g}\mathrm{e}})}^2\mathrm{d}t}{\int \mathrm{S}\mathrm{F}\mathrm{H}(t)\mathrm{d}t},$$(13)

representing the dispersion in the SFH of a clump. At all redshifts, the distribution peaks around log σ_SF/Myr = −1, with minimal spread, indicating that typical star formation episodes within clumps last approximately 0.1 Myr with little variability. However, each distribution shows a tail extending toward higher log σ_SF values, with some reaching up to 1.5. These higher values do not reflect prolonged star formation episodes but rather arise from secondary peaks in the SFHs. These secondary peaks have systematically low star formation rates, at least two orders of magnitude weaker than the primary peak, but with a similar duration. As an example, in systems with masses of 10⁴ M_⊙ the secondary peak corresponds to the formation of only ≃100 M_⊙, a negligible amount, comparable to that of a few massive stars. This small amount of mass is insufficient to be considered a substantial second star formation episode or accretion. Additionally, no clumps with masses exceeding 10⁴ M_⊙ exhibit multi-peaked SFH distributions, further underscoring that these secondary peaks represent minor contributions. This is a further clear indication that i) individual clumps rarely reaccrete gas to trigger subsequent star formation episodes and ii) that clumps are unlikely to result from mergers of multiple systems.

We conclude that mergers are not the primary driver of clump growth. This is supported by the star formation histories, which show isolated, short bursts, suggesting that clumps rarely re-accrete gas or undergo frequent mergers. It also reveals the way in which size–mass relations are typically populated. Thus, dissolution remains an important process affecting the low-mass end of the CMF.

Fig. 5 Scaling relations of clumps across redshifts. Top panels of each row: the CMF computed considering all (blue histogram) or only young (tage < 50 Myr, green histogram) clumps for the redshifts considered in this study. As a reference, the black dashed line indicates the $\frac{\mathrm{d}N}{\mathrm{d}{M}_{\star }}\propto M_{\star }^{-2}$ relation, corresponding to a scale free spectrum of structure formation. The dashed orange line shows, instead, the fit to the high mass end of the CMF (log M⋆ [M⊙] > 3.25), with the resulting slope indicated in the top right corner of each panel. Central panels of each row: Size-mass relation as a function of redshift (blue points with error bars). In each panel, the orange dashed line corresponds to a linear fit to the data points. The grey vertical lines and bands represent the region of space where tdiss ≤ 13 Gyr, with tdiss as in Equation (12). Each line corresponds to a different assumption on RG in Equation (12), specifically, from left to right, RG = 18 kpc, 12 kpc, 6 kpc, and 3 kpc, while VG = 220 km s−1 to mimic the tidal field the clumps would experience if orbiting around a Milky Way like galaxy (see Section 5.5 for detail). Bottom panels of each row: relation between the clumps stellar mass and Sérsic index. In the central and bottom panels, the red solid curve connects the median Re and n for each log M⋆ bin. Note that the panels at z = 6.14 cover a different mass range than the others to highlight the presence of a massive clump with M⋆ ≃ 2 × 106 M⊙.

Fig. 6 Top panel: SFH of the simulation at redshift 6.14, computed using two different binning schemes. The red line represents the SFH with bins of 0.2 Myr, while the black line uses bins of 10 Myr. Middle panel; star formation histories of a selection of clumps at z = 7. Each clump is color-coded differently, representing a range of variable masses and sizes. Bottom panel: cumulative stellar mass as a function of time. Note: the cumulative mass represents the total stellar mass formed throughout the simulation and does not account for mass loss due to massive stars. Therefore, the effective mass at redshift 6.14 is smaller (see Table 1). The lower panels illustrate the distribution of log σSF, which represents the spread in the SFHs of all clumps (see Equation (13)). This is presented for three distinct redshift values: z = 13, z = 9, and z = 7.

Fig. 7 Size-mass relation at a redshift of 7, providing a more detailed view compared to Fig. 5. Clumps are color-coded by age, with blue indicating younger clumps and red representing older clumps. Clumps have been further grouped into age bins: clumps younger than 0.16 Gyr, between 0.16 Gyr and 0.31 Gyr and older than 0.31 Gyr, with each bin containing an equal number of objects. The isodensity contours for the three age distributions are also shown as dashed blue (young), yellow (intermediate-age) and red (old) lines. The grey vertical lines and bands represent the region of space where tdiss ≤ 13 Gyr, with tdiss the same as in Equation (12). Each line corresponds to a different assumption on RG in Equation (12), specifically, from left to right, RG = 18 kpc, 12 kpc, 6 kpc, 3 kpc, together with VG = 220 km s−1 to mimic the tidal field of MW like galaxy (see text for details).

6 Discussion

In this section, we discuss our results in the context of the implemented stellar feedback model (Section 6.1), examining how it influences the formation and evolution of stellar clumps. We also compare our results with other simulation-based studies (Section 6.2), providing a framework to assess the impact of different simulation set-ups on the properties of high-redshift stellar systems.

6.1 The feedback model

Among the simulations presented in C24, a star formation model with 100% efficiency is necessary to create conditions suitable for producing massive and compact star clusters, underscoring the crucial impact of this parameter. Here, we provide a theoretical justification and describe the physical mechanism that facilitates the formation of these structures when ϵ_SFE = 1.

Recent studies (Dekel et al. 2023, but see also Renzini 2023, 2025) suggest that under the appropriate conditions, stellar winds, radiative, and supernova feedback may be insufficient to effectively heat or expel gas from star-forming regions, allowing star formation to proceed with remarkably high efficiency. This “feedback-free” starburst (FFB) scenario requires a lowmetallicity environment, typically at high redshift, and a gas density above a critical threshold (3 × 10³ cm⁻³). Once this density is reached, the free-fall and cooling timescales become comparable to the delay in feedback effectiveness, enabling a rapid accumulation and conversion of gas into stars.

Although differing in many respects, the FFB scenario and our star formation and feedback models share a number of similarities in their resulting phenomenology. In our simulations, star formation is triggered when gas cools below 2 × 10⁴ K and reaches a density threshold that is higher than 10⁴ cm⁻³ at the achieved resolution, exceeding the value required by the model from Dekel et al. (2023). At these densities, the free-fall time decreases to τ_FF < 0.5 Myr, while the cooling time similarly shrinks to comparable values. Once stars form, since radiative feedback is not accounted for, feedback occurs through two channels only: stellar winds and supernovae. Supernovae become effective only after about 6 Myr, corresponding to the lifetime of the most massive star (40 M_⊙) considered in the star formation model by C24. At this point, the differences in our model compared to the FFB framework are twofold: 1) stellar winds inject energy and mass continuously; 2) star formation efficiency remains fixed and does not depend on gas conditions.

For the FFB model to operate at this stage, a delayed onset of stellar winds would be necessary, which would automatically enhance the efficiency of star formation. In our case, setting ϵ_SFE = 1 accounts for this lack of self-consistency, enabling a rapid conversion of gas that proceeds until the injection of energy by winds is enough to adequately heat the remaining gas and bring the star formation episode to a halt. In our implementation, the thermal energy released by the winds increases rapidly as new stars are formed; because this energy is substantial even at low metallicity, it effectively heats the surrounding gas on a short timescale.

Figure 8 shows density and temperature maps for selected clumps at redshifts z = 15 (top rows) and z = 10.5 (bottom rows) after the star formation episode in the clump is quenched. Each slice is taken on a plane passing through the center of the target clump and young clumps were chosen to capture the feedback impact shortly after formation. Isodensity contours indicate projected stellar density along the corresponding line of sight. All clumps are embedded within large, high-temperature (>10⁷ K) and low-density (log n[cm⁻³] < −1) bubbles of gas recently heated by stellar winds and supernovae. In the two upper-left clumps, which are younger than 6 Myr, feedback exclusively comes from stellar winds, while in the other cases, both winds and supernovae contribute, as seen in the larger surrounding bubbles. This rapid, efficient star formation followed by intense heating drives the gas to a set of specific conditions (high temperature and low density) that is unfavorable for further star formation. Consequently, the star formation history (SFH) becomes highly episodic, with brief, intense bursts (Renzini 2023).

6.2 Comparison with previous works

This size–mass relation is a well-established trend, measured in observations (Adamo et al. 2024) and, to some extent, reproduced by simulations. Ma et al. (2020) analyze high-resolution cosmological simulations at redshifts greater than 5 to study the formation of dense and bound stellar clumps at high redshift, possibly progenitor of present-day GCs. In their work the authors focused on the evolution of halos that reach virial masses of 10¹⁰ M_⊙, 10¹ M_⊙, and 10¹² M_⊙ by a redshift of 5, considering different simulation resolutions as well. A key difference between our approaches is that Ma et al. (2020) did not model the formation of individual stars, but stellar particles instead (albeit of low mass) to represent stellar populations. However, as in our case, they adopted a high star formation efficiency, even though the specific implementation of their star formation differs slightly from ours. Ma et al. (2020) find that all their simulated halos are able to form bound star clusters, which then follow a similar size–mass relation. However, their results demonstrate a great sensitivity to resolution, which causes the minimum mass and the typical sizes of the clusters to be resolution-dependent. Only their highest resolution simulation is comparable with ours and, in that case, for clumps with masses in the range of log M_⋆ [M_⊙] = 3.5–4.5, they reported sizes between 2 pc and 100 pc (as shown in their Fig. 13), which overlap with the sizes we measure for clumps within the same mass range. However, their size distribution shows a significantly larger scatter compared to ours. In spite of this, both studies are in agreement with respect to the CMF, reporting a power-law slope close to −2.

In another study, by means of hydrodynamical, cosmological simulations at subparsec resolution, Garcia et al. (2023) investigate the formation and evolution of star clusters in galaxies in the typical mass regime of dwarfs (M_vir ≃ 2 × 10⁸ at z = 8). These simulations include radiative feedback and do not model individual stars as we do, while exploring the effects of adopting two distinct star formation efficiencies: 0.35 (low) and 0.7 (high). The stellar clusters formed in their simulations span a lower mass range compared to ours, from approximately a few 10² M_⊙ to 10⁴ M_⊙. As discussed by the authors, the CMF in their simulations changes significantly in response to repeated episodes of star formation, with the CMF flattening after each intense burst of new star formation. In their low-SFE run, Garcia et al. (2023) found a CMF slope similar to ours, but in their high-SFE case, the slope becomes shallower: around −0.6. This presents a key divergence from our results, as in our simulations, comparable to their high-SFE run, the CMF retains a steep slope of approximately −2. In their high-SFE run, the authors discussed how the bound clusters formed are systematically denser compared to the low-SFE case. This aligns with C24, where bound star cluster densities reach up to almost 10⁴ M_⊙ pc⁻² only in the high ϵ_SFE run, while they are much less dense and more diffuse in low ϵ_SFE simulations. As for the size–mass relation, Garcia et al. (2023) also observed a clear positive correlation. Additionally, their work finds a correlation between the size–mass relation and the ages of the clumps, indicating that younger clumps tend to populate the region of higher mass and larger size. This age correlation is consistent with our findings, further supporting the relationship between clump size, mass, and age.

Gutcke (2024) presented high-resolution cosmological simulations aimed at exploring scenarios for GC formation. The simulations focus on dwarf galaxies with stellar masses, at z = 0, ranging from 10⁶ M_⊙ to 10⁷ M_⊙, closely matching the properties of Local Group dwarf galaxies. The study reveals that clusters host ancient stellar populations and exhibit short, episodic star formation histories, occasionally marked by the presence of multiple stellar generations. Our size–mass relation is consistent with that of the author, showing a very similar correlation; although we note that Gutcke (2024) forms clusters with a generally smaller average mass than ours. This may be attributed to the lower mass of the dark matter halo in their simulation. This author also found very small fractions of dark matter relative to stellar mass in the stellar clusters formed, accompanied by a discussion of how dark matter is lost through tidal stripping. While also the majority of our clusters live in environments free from dark matter, in our case, this absence is attributable to a formation without dark matter rather than a loss over time. In addition, while the typical SFHs in Gutcke (2024) and our clumps are characterized by short bursts of star formation, our clumps do not show the capability to re-accrete gas and they also lack multimodal distributions, which diverges from the findings of Gutcke (2024).

Fig. 8 Top panels: slice of gas density passing through the center of four selected clumps (first row). Darker colors indicate regions of higher density, while lighter colors correspond to lower density areas. The second row displays the corresponding temperature maps. Superimposed on each map are curves of constant projected stellar density, shown as white lines (top) and black lines (bottom). All clumps in the first row are observed at redshift z = 15.5. Bottom panels: same as the top panels, but for four clumps at z = 10.5. In the top right part of each panel, we have indicated the age of the clump.

7 Summary and conclusions

In this study, we explore the formation of compact, high-redshift star-forming clumps using high-resolution cosmological zoom-in simulations. This simulation is part of the SIEGE project (Calura et al. 2022), a suite designed to investigate the conditions of the formation of compact star clusters that are possible progenitors of today’s GCs. Specifically, this simulation is drawn from the larger suite described by C24, which examines the effects of feedback models, star formation efficiencies, and IMF variations on the structural and dynamical properties of young star clusters in high-redshift galaxies. The simulation targeted a dwarf galaxy halo with a virial mass of 4 × 10⁹ M_⊙ at z = 6.14, achieving subparsec spatial resolution. This exceptional resolution, combined with individual star formation sampling from the IMF and an enforced 100% star formation efficiency, has enabled us to investigate the conditions and mechanisms that govern star formation in the extreme environments characteristic of high-redshift protogalaxies. In this work, we specifically investigate the properties of the stellar cluster population formed in our simulation, focusing on the evolution of their characteristics and scaling relations with redshift.

The high-redshift environment and the star formation and feedback model probed by our simulation naturally favor the formation of bound stellar clumps, with masses spanning the range 10³ M_⊙ to 10⁶ M_⊙ and with characteristic sizes on the order of 1–3 parsecs, along with surface densities reaching up to almost 10⁴ M_⊙ pc⁻². The spatial resolution of the simulation (≃0.2–0.5 pc across the redshift range) is comparable to the size of the smallest stellar clumps formed. As such, the clump sizes reported should be regarded as upper limits. While higher-resolution simulations could refine the structural characterization of clumps and potentially reveal more compact ones, such simulations are currently computationally unfeasible. Conversely, lower-resolution runs would artificially inflate clump sizes due to degraded gravitational accuracy, thereby offering no meaningful convergence test. Remarkably, nearly all clumps are found to have minimal dark matter, with a ratio of dark-to-stellar mass within three effective radii below 1 in most cases. This finding suggests that such dense, compact structures can originate independently of dense dark matter halos – initially forming as bound star clusters, rather than transient associations. Their eventual dissolution occurs on longer timescales, primarily due to interactions with the host galaxy. We reproduced the typical behavior of the CMF, which follows a power law of slope approximately −2. All clumps adhere to a mass-size relation that exhibits little evolution with redshift: smaller clumps tend to be less massive than larger clumps. The size mass relation is in agreement with results from other works. We found a clear trend in the size–mass relation with clumps age: younger clumps generally populate the higher mass and size regions. Clumps at higher redshift appear more compact, showing an evolution in surface density as we move to earlier cosmic times. The SFH of these clumps is consistently mono-modal, indicating well-defined ages and formation in single, rapid star formation episodes of a typical duration of ≃0.1 Myr. We observe no multimodal SFHs, implying a limited capacity for gas re-accretion and subsequent star formation bursts in these environments.

Our ability to resolve subparsec scales and incorporate feedback from individual stars is essential for capturing the small-scale processes governing clump formation and evolution in high-redshift environments. This study provides a foundation for understanding the physical conditions under which GCs evolve, while demonstrating the critical role of star formation efficiency as a key parameter in formation models of high-redshift clumps.

Acknowledgements

This paper is supported by the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations). The research activities described in this paper have been co-funded by the European Union – NextGeneration EU within PRIN 2022 project no. 20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch. We acknowledge support from the INAF Minigrant ‘Clumps at cosmological distance: revealing their formation, nature, and evolution’ (Ob. Fu. 1.05.23.04.01). We acknowledge PRACE for awarding us access to Discoverer at Sofia Tech Park, Bulgaria. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. We acknowledge EuroHPC JU for awarding the project IDs EHPC-REG-2021R0052 and EHPC-REG-2024R01-042 access to DISCOVERER at the Sofia Tech Park, Bulgaria. AL acknowledges support from PRIN MUR “2022935STW” funded by European Union-Next Generation EU, Missione 4 Componente 2 CUP C53D23000950006.

References

All Tables

All Figures

Fig. 1 Evolution with redshift of the projected dark matter density within the high-resolution simulation domain. Top, middle, and bottom panels show the projected dark matter density distribution in the central region of the simulation box at a redshift of z = 10.5, 8, and 7, respectively. In all panels, brighter colors indicate higher-density regions. The stellar projected density maps corresponding to the same regions are overlaid for comparison, allowing for a visual examination of the relationship between dark matter and stellar structures across different cosmic epochs. The stellar clumps formed have been highlighted with black circles.

In the text

Fig. 2 Examples of stellar density distributions of star-forming regions taken at three different redshifts. The top-right, top-left, and lower maps correspond to z = 10.5,7, and 6.14, respectively. The star-forming systems shown in the top-right and top-left panels correspond to the stellar counterparts of the most massive halos visible in the top and bottom panels of Fig. 1. In each panel, the small insets display the star formation histories of the corresponding regions. All clumps identified within the time windows color-coded in the small insets are marked with circles and displayed with the same colors on the corresponding projected stellar density map.

In the text

Fig. 3 Clump projected density distributions. Top panels: projected stellar density profiles (black squares with errobars) for two clumps selected at z = 15.5. Each set of three panels shows the projected density along different lines of sight: along the z, y, and x axes from left to right, respectively (the x z plane is the plane of Fig. 1). The blue curves represent the best-fitting Sérsic + background models, while the vertical lines indicate the position of the effective radius, Re. In the middle and right panels of each triplet we added (with a grey line) the Sérsic + background model shown in the left panel to highlight differences among projections. The curves and the binned profiles shown are those obtained at the end of the two step fitting procedure (see Section 4). The small insets display the projected density map along the same direction as in the main panel, with the blue circle centered at center of the clump and having a radius equal to Re. The clump in the left panel has (M⋆, Re, n) = (1.25 × 104 M⊙, 0.67 pc, 1.24), while the one in the right panel has (M⋆, Re, n) = (1.79 × 104 M⊙, 0.76 pc, 1.47). Middle panels: same as the top panels, but for clumps at z = 10.5. The clumps have (M⋆, Re, n) = (2.88 × 103 M⊙, 1.00 pc, 1.63) – left – and (M⋆, Re, n) = (3.04 × 104, 2.13 pc, 2.69) – right. Bottom panels: Same as the top and middle panels, but for clumps at z = 8, with (M⋆, Re, n) = (3.03 × 103 M⊙, 1.57 pc, 1.62), and (M⋆, Re, n) = (1.83 × 104 M⊙, 1.71 pc, 2.27) on the right. In each panel, the orange shaded area shows the simulation resolution (Δx) at that redshift, Δx ≡ 5h−1/[2l(1 + z)], with l = 21 and h = 0.703.

In the text

Fig. 4 Clump parameters across redshifts. The triplets of panels illustrate: the relation between the dark-matter-to-stellar-mass ratio and clump ages (top-left), the distributions of the dark matter mass within three effective radii relative to the total stellar mass, Mdm (3Re)/M⋆ (top-right), the distributions of the boundedness parameter Π for clumps younger than 50 Myr (middle-left), the distributions of the effective radius Re (middle right), the distributions of the Sérsic index n (bottom left), and the distributions of the average surface density ${\mathrm{\Sigma }}_{\mathrm{e}}\equiv \frac{M_{\star }}{2\pi R_{\mathrm{e}}^2}$ (bottom-right). Each panel within a triplet represents distributions at different redshifts, with z = 13 shown in blue (top), z = 9 in purple (middle), and z = 7 in orange (bottom). The colored bands in the top left triplet of panels highlight clumps younger than 50 Myr. These clumps are used to build the log Π distributions in the middle-left triplet of panels. For clarity, only three of the eight redshifts analyzed in this work are displayed. The vertical dashed line in all triplets (except for the top left) show the median value of the corresponding distribution.

In the text

Fig. 5 Scaling relations of clumps across redshifts. Top panels of each row: the CMF computed considering all (blue histogram) or only young (tage < 50 Myr, green histogram) clumps for the redshifts considered in this study. As a reference, the black dashed line indicates the $\frac{\mathrm{d}N}{\mathrm{d}{M}_{\star }}\propto M_{\star }^{-2}$ relation, corresponding to a scale free spectrum of structure formation. The dashed orange line shows, instead, the fit to the high mass end of the CMF (log M⋆ [M⊙] > 3.25), with the resulting slope indicated in the top right corner of each panel. Central panels of each row: Size-mass relation as a function of redshift (blue points with error bars). In each panel, the orange dashed line corresponds to a linear fit to the data points. The grey vertical lines and bands represent the region of space where tdiss ≤ 13 Gyr, with tdiss as in Equation (12). Each line corresponds to a different assumption on RG in Equation (12), specifically, from left to right, RG = 18 kpc, 12 kpc, 6 kpc, and 3 kpc, while VG = 220 km s−1 to mimic the tidal field the clumps would experience if orbiting around a Milky Way like galaxy (see Section 5.5 for detail). Bottom panels of each row: relation between the clumps stellar mass and Sérsic index. In the central and bottom panels, the red solid curve connects the median Re and n for each log M⋆ bin. Note that the panels at z = 6.14 cover a different mass range than the others to highlight the presence of a massive clump with M⋆ ≃ 2 × 106 M⊙.

In the text

Fig. 6 Top panel: SFH of the simulation at redshift 6.14, computed using two different binning schemes. The red line represents the SFH with bins of 0.2 Myr, while the black line uses bins of 10 Myr. Middle panel; star formation histories of a selection of clumps at z = 7. Each clump is color-coded differently, representing a range of variable masses and sizes. Bottom panel: cumulative stellar mass as a function of time. Note: the cumulative mass represents the total stellar mass formed throughout the simulation and does not account for mass loss due to massive stars. Therefore, the effective mass at redshift 6.14 is smaller (see Table 1). The lower panels illustrate the distribution of log σSF, which represents the spread in the SFHs of all clumps (see Equation (13)). This is presented for three distinct redshift values: z = 13, z = 9, and z = 7.

In the text

Fig. 7 Size-mass relation at a redshift of 7, providing a more detailed view compared to Fig. 5. Clumps are color-coded by age, with blue indicating younger clumps and red representing older clumps. Clumps have been further grouped into age bins: clumps younger than 0.16 Gyr, between 0.16 Gyr and 0.31 Gyr and older than 0.31 Gyr, with each bin containing an equal number of objects. The isodensity contours for the three age distributions are also shown as dashed blue (young), yellow (intermediate-age) and red (old) lines. The grey vertical lines and bands represent the region of space where tdiss ≤ 13 Gyr, with tdiss the same as in Equation (12). Each line corresponds to a different assumption on RG in Equation (12), specifically, from left to right, RG = 18 kpc, 12 kpc, 6 kpc, 3 kpc, together with VG = 220 km s−1 to mimic the tidal field of MW like galaxy (see text for details).

In the text

Fig. 8 Top panels: slice of gas density passing through the center of four selected clumps (first row). Darker colors indicate regions of higher density, while lighter colors correspond to lower density areas. The second row displays the corresponding temperature maps. Superimposed on each map are curves of constant projected stellar density, shown as white lines (top) and black lines (bottom). All clumps in the first row are observed at redshift z = 15.5. Bottom panels: same as the top panels, but for four clumps at z = 10.5. In the top right part of each panel, we have indicated the age of the clump.

In the text