© The Authors 2026

1 Introduction

Whether stars in stellar systems form through stochastic sampling or through ordered, physically regulated processes has direct consequences for how stellar populations, cluster demographics, and feedback are interpreted. Most analytic models and population synthesis approaches implicitly assume that stellar masses are drawn randomly from a universal initial mass function (IMF; Kroupa 2001b; Chabrier 2003). However, observational and theoretical work challenges this view. Weidner et al. (2013) and Yan et al. (2023) show that the observed relation between the stellar mass of the embedded star cluster, M_ecl, and the maximum stellar mass, m_max, is inconsistent with purely random sampling, indicating physical regulation during cluster formation. The relation between m_max and the mass of the molecular cloud clump was noticed by Larson (1982) and formulated as a power-law relation in Larson (2003). The existence of a well-defined relation between m_max and M_ecl for clusters born in molecular clumps was quantified for the first time by Weidner & Kroupa (2006). They showed with multiple Monte Carlo experiments that observational data of young clusters are in better agreement with an ordered formation than a random one. This has been interpreted as evidence of deterministic IMF sampling, of which optimal sampling is a specific case (Kroupa et al. 2013; Schulz et al. 2015).

Optimal sampling is a specific deterministic scheme in which stellar masses are assigned so that the IMF is matched exactly in all mass intervals. Deterministic sampling refers to any sampling scheme in which a given IMF and total stellar content uniquely and reproducibly determine the resulting set of stellar masses, implying a highly regulated form of star formation with no intrinsic sampling noise. Random sampling treats the IMF as a probability distribution from which stellar masses are drawn independently, implying that star formation is intrinsically stochastic and that finite stellar populations exhibit unavoidable statistical fluctuations about the IMF (Kroupa 2008). Weidner & Kroupa (2006) ruled out random sampling with 99.9% confidence, which was confirmed in Weidner et al. (2013). Using observational data of the young star-forming regions Taurus, Lupus3, ChaI, and IC348 (Kirk & Myers 2012), young stellar clusters in the Large Magellanic Cloud (Stephens et al. 2017), the star-forming regions Taurus-Auriga and L1641 in Orion (Weidner et al. 2013), and new observational data from the VVV survey (Ramírez Alegría et al. 2016), Yan et al. (2017) show that the m_max–M_ecl relation is consistent with deterministic sampling approaches. Random sampling can also reproduce the overall trend, but typically with a larger intrinsic scatter than observed (Yan et al. 2023). Furthermore, deterministic and optimal sampling are consistent with a highly self-regulated picture of star formation (Yan et al. 2017; Vázquez-Semadeni et al. 2024).

The sampling method also affects our understanding of galaxies whose stellar content is the sum of the stars in the galaxy’s embedded star clusters, in which the majority of stars form (Kroupa & Weidner 2003). This concept is included in the theory of the integrated galaxy-wide IMF, which comprises the distribution of stellar remnants, the number of supernovae, and the chemical enrichment history of a galaxy (Weidner & Kroupa 2005; Jeřábková et al. 2018; Haslbauer et al. 2024; Zonoozi et al. 2025). This theory also has several cosmological implications, which are listed in Kroupa et al. (2026), and it solves multiple problems, such as the short star formation timescales of elliptical galaxies, named the downsizing problem (Yan et al. 2021), with significant implications for our understanding of the cosmic microwave background (Gjergo & Kroupa 2025).

Despite evidence of the m_max–M_ecl relation across multiple associations and spatial scales, several studies have challenged it (e.g. Andrews et al. 2013; Jung et al. 2023). Weidner et al. (2013) argued that the reported contradictions can arise from the interpretation of the observational data given observational biases. By revisiting the analysis of Andrews et al. (2013), Weidner et al. (2014) defended the m_max–M_ecl relation. The present work revisits the interpretation presented by Jung et al. (2023, hereafter J23). In their analysis, the authors report no clear evidence for the m_max–M_ecl relation in their observed megaparsec-distant galaxies of the Legacy ExtraGalactic UV Survey (LEGUS). Young, massive stars emit high-energy photons that ionise the surrounding gas. Hydrogen recombination produces a prominent Hα emission line, which is therefore a useful tracer of strongly ionising stellar populations (Pflamm-Altenburg et al. 2007). Assuming a canonical IMF and a maximum stellar mass of 120 M_⊙, the approach by J23 was to analyse the dependence of the upper end of the IMF on the cluster mass by normalising the measured Hα luminosity, L_Hα, by M_ecl and plotting it against M_ecl¹. The result is displayed in their Fig. 4. In that diagnostic plot, the clusters appear broadly consistent with a randomly populated distribution, from which J23 conclude that the m_max–M_ecl relation is not evident in their observed galaxies.

The goal of this work is to test which cluster-sampling process is most consistent with observations. We optimally and randomly sampled star clusters and determined their Hα luminosities for both sampling cases. We then applied evolutionary processes (stellar mass evolution and dynamical ejections) and considered observational uncertainties (age uncertainty and contamination by field stars). Finally, we compared the resulting cluster distributions to the observations.

The concepts that are relevant for this work are described in Sect. 2. The results are shown in Sect. 3 and are discussed in Sect. 4. This work is summarised in Sect. 5.

2 Methods

In this section we describe the framework used to model the stellar content and the resulting Hα luminosity of young star clusters. The primary goal is to predict the ionising output of a cluster of given mass and age by explicitly accounting for the stochastic and systematic effects associated with sampling the stellar IMF and subsequent evolutionary processes. The model combines a variable IMF prescription, different IMF sampling schemes, stellar evolution–based mass and luminosity estimates, and a treatment of observational and physical uncertainties. Together, these ingredients allow us to quantify how cluster mass, age, and IMF sampling affect the inferred Hα emission and the observable stellar population. Twenty embedded cluster masses, M_ecl, were chosen uniformly distributed over log-mass in the range 10^2.5 M_⊙ to 10^4.5 M_⊙. We note here that we merely needed to cover the approximate range of clusters covered observationally by J23 but with a larger mass range to account for clusters that move out of the observational range due to evolutionary processes and the biases studied here. The ensemble of embedded cluster masses thus comprises the following values in M_⊙: 316, 428, 579, 784, 1062, 1438, 1947, 2636, 3569, 4832, 6543, 8858, 11993, 16237, 21983, 29763, 40296, 54555, 73861, and 100000. For the final comparison to J23, we restricted the model to the mass range covered by the observed sample and down-sampled the highest-mass model clusters to reflect the fact that J23 included only a few clusters at the high-mass end.

2.1 The IMF

The initial stellar mass function was introduced as a one-part power law, dN = ξ(m) dm = km^−α dm, where dN is the number of stars in the mass interval between m and m + dm, k is a normalisation constant, and α is the slope of the mass function. Salpeter (1955) derived the IMF in the mass range 0.4 ≤ m/M_⊙ ≤ 10 with a slope of α = 2.35 from the luminosity function of main-sequence stars in the solar neighbourhood.

Based on an extensive star-count analysis of the Galactic field by Kroupa et al. (1993) and observations of young star clusters and OB associations in the Milky Way and the Magellanic Clouds, Kroupa (2001a) concluded that the IMF is a multi-part power law; the data support a three-part power law (the canonical IMF): $$\xi (m)=k\{\begin{array}{ll}{m}^{-{\alpha }_1},& m_{min}<m\le m_1,\\ k_2\cdot m^{-{\alpha }_2},& m_1<m\le m_2,\\ k_3\cdot m^{-{\alpha }_3},& m_2<m\le m_{max},\end{array}$$(1)

where k_i has to be determined such that the IMF fulfils continuity and k is needed for normalisation. Furthermore, m_min = 0.08 M_⊙ and m_max are the lowest and maximum stellar mass, respectively. Studying the rich star clusters R136 and Arches, masses up to 150 M_⊙ could be observed (Weidner & Kroupa 2004; Figer 2005), which will be adopted in this work as the maximum physical stellar mass. More massive stars are likely mergers (Banerjee et al. 2012) and brown dwarfs appear as objects associated with the formation of stars (Kroupa et al. 2026 for further information). In the canonical IMF, α₁ = 1.3 and α₂ = α₃ = 2.3.

The systematic change of the IMF is observable in extreme star-forming environments and metallicities different compared to the Milky Way’s as shown in Dabringhausen et al. (2008, 2009, 2012). Marks et al. (2012) found a ’fundamental plane’ that describes the change of the IMF with density and metallicity. The authors also find that the density of the molecular cloud clump has a bigger impact on the upper end of the IMF than the metallicity.

In Yan et al. (2017), Eqs. (2) and (3) were empirically determined and Δα = 63 is suggested for Z_⊙ = 0.0142 (Yan et al. 2021), $${\alpha }_1=0.3+\mathrm{\Delta }\alpha \cdot (\mathrm{Z}-{\mathrm{Z}}_{\odot })0.08\le m/M_{\odot }<0.5,$$(2) $${\alpha }_2=1.3+\mathrm{\Delta }\alpha \cdot (\mathrm{Z}-{\mathrm{Z}}_{\odot })0.5\le m/M_{\odot }<\mathrm{1.0.}$$(3)

The dependence of the power-law index for high-mass stars was determined empirically by Marks et al. (2012). The formulation by Yan et al. (2017) that is used in this work is $${\alpha }_3=\{\begin{array}{ll}2.3,& x<-0.87\\ -0.41x+1.94,& x>-0.87\end{array}1\le \frac{m}{M_{\odot }}<m_{max},$$(4)

with x = −0.14[Z] + 0.99 log₁₀(ρ_cl/(10⁶M_⊙pc⁻³)), where ρ_cl is the clump density (final stars plus residual gas assuming a star formation efficiency of 1/3) and [Z] = log₁₀(Z/Z_⊙) ≈ [Z/X] with the initial metal-mass fraction, Z, and hydrogen-mass fraction, X, in stars.

2.2 Sampling of star clusters

To obtain a finite stellar population, the IMF must be discretised. Here we considered two widely used discretisation schemes: random sampling (Sect. 2.2.1) and optimal sampling (Sect. 2.2.2). In optimally sampled embedded clusters, stellar masses are assigned such that the discretised population matches the IMF exactly in each mass interval; by construction this minimises sampling noise and avoids gaps (Kroupa et al. 2013; Yan et al. 2017). In contrast, randomly sampled clusters treat the IMF as a probability distribution and therefore show finite-population fluctuations (including gaps), and can occasionally contain more massive stars at fixed M_ecl compared to the optimal-sampling case (see Fig. 1 of Kroupa et al. 2013).

2.2.1 Random sampling

For random sampling, the IMF is treated as a probability density function, p(m) ∝ ξ(m). Using a generating function, stellar masses can be drawn randomly from the IMF as described in Kroupa (2008). First, the normalisation constant needs to be determined using the following condition: $$1={\int }_{m_{min}}^{m_{max}}p(m)\mathrm{d}m.$$

From the continuity condition $m_1^{-{\alpha }_1}=k_2{m}_2^{-{\alpha }_2}$ and $k_2{m}_2^{-{\alpha }_2}=k_3{m}_3^{-{\alpha }_3},k_2$ and k₃ can be determined.

After defining the generating functions with the IMF slopes in Eqs. (2)–(4) and a random number between 0 and 1, the corresponding stellar mass was obtained as described in Kroupa (2008). Based on this concept, we developed a Python code² that draws stellar masses until the target M_ecl is reached, and then removes the final draw to avoid exceeding the target mass. Consequently, the realised stellar mass is slightly below the initial target M_ecl; throughout the analysis we therefore used the realised cluster mass, computed by summing the drawn stellar masses.

2.2.2 Optimal sampling

In optimal sampling, the IMF is assumed to be an optimal distribution function where the number of stars in any stellar mass range matches exactly the integration of the given IMF in that mass range (Kroupa et al. 2013). Codes for optimal sampling are provided by Yan et al. (2017)³ and Gjergo et al. (2026)⁴, based on the mathematical formulation of optimal sampling by Schulz et al. (2015).

To optimally sample the stellar IMF in an embedded cluster, the normalisation constant and m_max need to be determined with the following conditions: the mass conservation law $M_{\text{ecl}}={\int }_{m_{\text{min}}}^{m_{\text{max}}}m\xi (m)\mathrm{d}m$ and the optimal-sampling normalisation condition $1={\int }_{m_{\text{max}}}^{150M_{\odot }}\xi (m)\mathrm{d}m$. The different integral boundaries, m_i, are determined by summing up all stars as described in Eq. (5), where each integral must be one and N_tot is the calculated total number of stars, $$\begin{array}{rl}{N}_{\mathrm{t}\mathrm{o}\mathrm{t}}& ={\int }_{m_{min}}^{m_{max}}\xi (m)\mathrm{d}m\\ & ={\int }_{m_2}^{m_{max}}\xi (m)\mathrm{d}m+\dots +{\int }_{m_{N_{\mathrm{t}\mathrm{o}\mathrm{t}}}}^{m_{N_{\mathrm{t}\mathrm{o}\mathrm{t}}-1}}\xi (m)\mathrm{d}m+{\int }_{m_{min}}^{m_{N_{\mathrm{t}\mathrm{o}\mathrm{t}}}}\xi (m)\mathrm{d}m.\end{array}$$(5)

The calculation starts with the most massive star. So with ${\int }_{m_2}^{m_{\text{max\hspace{0.17em}}}}\xi (m)\mathrm{d}m=1$ and the before deduced m_max, m₂ can be calculated. Then, the remaining m_i are determined iteratively. With $m_i={\int }_{m_{i+1}}^{m_i}m\xi (m)\mathrm{d}m$, the corresponding mass of the star is calculated.

2.3 Stellar evolution

We assumed that all stars in an embedded cluster have the same initial age. We accounted for stellar mass evolution with time using evolutionary tracks⁵. They were selected to start from the pre-main sequence and continue until either the first thermal pulsing or C-ignition from PARSEC (Bressan et al. 2012) version 1.2 S, which is available for 0.0001 ≤ Z ≤ 0.02 in the mass range 0.1 ≤ m/M_⊙ < 350, for 0.03 ≤ Z ≤ 0.04 in the mass range 0.1 ≤ m/M_⊙ < 150, and for Z = 0.06 in the mass range 0.1 ≤ m/M_⊙ < 20 (cf. Tang et al. 2014 for 0.001 ≤ Z ≤ 0.004 and Chen et al. 2015 for other Z), with revised and calibrated surface boundary conditions in low-mass dwarfs provided in Chen et al. (2014). Furthermore, the thermal pulsing-asymptotic giant branch evolution from the first thermal pulse to the total loss of the envelope is added with COLIBRI S_37 (Pastorelli et al. 2020) tracks for 0.008 ≤ Z ≤ 0.02, COLIBRI S_35 tracks (Pastorelli et al. 2019) for 0.0005 ≤ Z ≤ 0.006 and COLIBRI PR16 (Marigo et al. 2017) for Z ≤ 0.0002 and Z ≥ 0.03.

The tracks were downloaded for ages between 1 and 10 Myr and Z = {0.004, 0.008, Z_⊙ ≈ 0.014}, where 0.014 is the solar metallicity according to Asplund et al. (2009). With interpolation, the mass of a star at a given age, m, is connected to the initial mass, m_ini.

2.4 Ho luminosity

The Hα luminosity was computed using the galaxy evolution code Pégase (Fioc & Rocca-Volmerange 1999; Fioc & Rocca-Volmerange 2019), which models the spectra and chemical evolution of galaxies across a wide wavelength range⁶. We computed synthetic spectra by setting up simple stellar populations (SSPs), where all stars have the same mass and age. The technical details are described in Appendix A. The obtained Hα luminosity of the SSP was normalised to the one of a single star. To do so, the number of objects in the SSPs, N_tot, needed to be obtained with the given IMF and the given stellar mass, which was normalised to 1 M_⊙ by Pégase.

2.5 Dynamical ejections

For clusters younger than 4 Myr, gas expulsion has not significantly changed M_ecl (Weidner et al. 2010), but the dynamical ejection of stars may have an impact. Depending on the initial radius, clusters can eject their most massive stars within the first 3 Myr (Oh & Kroupa 2012). To assess the influence of dynamical ejections on the m_max–M_ecl relation, we adopted the results of Oh et al. (2015) and Oh & Kroupa (2012, 2016).

Oh & Kroupa (2012) performed NBODY6 (Aarseth 1999) simulations of young star clusters with masses in the range M_ecl = 10–10^3.5 M_⊙, half-mass radii of r_h = {0.3, 0.8} pc, and binary fractions of 0 (all single stars) and 1 (all binaries). The binaries were randomly paired or ordered paired; the latter means that stars heavier than 5 M_⊙ were paired to binaries with mass ratios close to unity. Lower-mass stars were randomly paired (for a discussion, see Kroupa et al. 2026 and Kroupa 2025). In Oh & Kroupa (2012), and the positions and velocities of stars in the clusters were generated according to the Plummer model, which is explained in Aarseth et al. (1974). For stellar evolution, the stellar evolution library (Hurley et al. 2000) in the NBODY6 code was used. Oh & Kroupa (2016) analysed the simulations performed in Oh & Kroupa (2012) and Oh et al. (2015) for M_ecl = 10⁴ M_⊙ and 10^4.5 M_⊙.

Referring to Goodwin & Kroupa (2005), most stars form in systems of two or three stars. Furthermore, Sana et al. (2008) studied O-type stars in NGC 6231, finding that random pairing is unlikely for O-type stars and that O+OB binaries are favoured. This is confirmed in Sana et al. (2009) and other publications. Furthermore, star clusters are commonly initially mass-segregated (Bonnell & Davies 1998; Gouliermis et al. 2004; Chen et al. 2007; Pavlík et al. 2019), which is accounted for in the simulations by Oh & Kroupa (2012) and Oh et al. (2015).

Marks & Kroupa (2012) performed an inverse dynamical population synthesis. By studying binary populations together with constraints on the birth densities of globular clusters obtained in Marks & Kroupa (2010), a correlation between M_ecl and the half-mass radius r_h of the form $\frac{r_{\mathrm{h}}}{\mathrm{p}\mathrm{c}}={0.10}_{-0.04}^{+0.07}$. ${\left(\frac{M_{\text{ecl\hspace{0.17em}}}}{M_{\odot }}\right)}^{0.13\pm 0.04}$ was found. With this, r_h of the clusters considered in this analysis is estimated to be 0.3 pc within 1σ confidence. We therefore adopted model MS30P of Oh et al. (2015), which assumes initial mass segregation, a half-mass radius of r_h ≈ 0.3 pc, and ordered pairing for massive binaries.

The ejection fraction depends on the spectral type and, hence, stellar mass (Oh & Kroupa 2016, Fig. 1). The division of the spectral types is displayed in Table B.1. The mean ejection fraction of O-star systems, i.e. of multiple stars, at an age of $3\mathrm{M}\mathrm{y}\mathrm{r},⟨f_{\mathrm{e}\mathrm{j},\mathrm{O}-\text{sys}}⟩=\frac{N_{\mathrm{e}\mathrm{j},\mathrm{O}-\text{sys}}}{N_{\mathrm{O}-\text{sys}}}$, with the number of ejected O-star systems, N_ej,O–sys, and the number of O-star systems in the cluster, N_O–sys, is listed in Table B.2 according to Oh et al. (2015, their Table A1). We note that the notation used in this work differs slightly from that in the original reference. A star or star system is classified as ejected if its distance from the cluster centre is greater than 3 × r_h of the cluster at 3 Myr, and its velocity is larger than the escape velocity at that radius. With a Gaussian function of the form $f_{\mathrm{e}\mathrm{j},\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}\approx 0.23\cdot \exp (-\frac{1}{2}{\left(\frac{\log \left(M_{\mathrm{e}\mathrm{c}\mathrm{l}}/M_{\odot }\right)-3.5}{-0.70}\right)}^2)$, the ejection fractions of O-star systems for different cluster masses were inferred. Furthermore, the uncertainties were taken into consideration by adding the values of the standard deviation, which are listed in Table B.2 and were determined from the values of all data points, and calculating the standard deviation. In the case of multiple points with the same value where the determination of the number of points is difficult, the number is estimated to yield the number of data points (N_run in Table B.2) and the mean value to be in agreement with ⟨f_ej,O–sys⟩ in Table B.2. In order to determine σ_{⟨fej,O–sys⟩} for a M_ecl that is not listed in Table B.2, the value of the closest M_ecl for which a value for σ_{⟨fej,O–sys⟩} is listed is assumed. Thus, the ejection fraction can be increased within σ confidence regions. In the case of a 2σ higher ejection fraction, σ_{⟨fej,O–sys⟩} is assumed to be twice as high.

Oh et al. (2015) present the averaged binary fraction of ejected O-star systems in their Fig. 11. The values are listed in the third column of our Table B.2. To obtain the binary fraction, we fitted a function of the form ⟨f_b,ej,O⟩ ≈ 32.07 · (M_ecl/M_⊙)^−0.65 to the data points of model MS30P. They found that the mass-ratio distribution does not significantly influence the ejection of massive stars as long as massive stars are paired with other massive stars in binaries. We therefore considered stars from the same spectral type for the ejection of systems but did not apply ordered pairing to the collection of ejected stars.

Assuming a constant ejection rate during ages from 1–3 Myr, the fraction of ejections at a specific age t can be determined with $$f_{\mathrm{t}}=\{\begin{array}{ll}(t/\mathrm{M}\mathrm{y}\mathrm{r})/3& \text{if}t\le 3\mathrm{M}\mathrm{y}\mathrm{r},\\ 1& \text{if}t>3\mathrm{M}\mathrm{y}\mathrm{r}.\end{array}$$

The number of ejected O stars is thus N_ej,O = N_ej,O–sys · (1 + f_b,ej,O) · f_t. Note that $f_{\mathrm{e}\mathrm{j},\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}=\frac{N_{\mathrm{e}\mathrm{j},\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}}{N_{\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}}(\mathrm{not}\frac{N_{\mathrm{e}\mathrm{j},\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}}{N_{\mathrm{o}}})$. If all stars are in binaries, then N_O–sys can be approximated as $N_{\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}\approx \frac{1}{2}{N}_{\mathrm{O}}$. So, the final ejection fraction of O-type stars is calculated as $$f_{\mathrm{e}\mathrm{j},\mathrm{O}}=\frac{N_{\mathrm{e}\mathrm{j},\mathrm{O}}}{N_{\mathrm{O}}}=\frac{1}{2}\cdot f_{\mathrm{e}\mathrm{j},\mathrm{O}-\mathrm{s}\mathrm{y}\mathrm{s}}\cdot (1+f_{\mathrm{b},\mathrm{e}\mathrm{j},\mathrm{O}})\cdot f_{\mathrm{t}}.$$

With the fraction conversion f_conv displayed in Table B.1, the ejection fractions of the other spectral types are determined. The values of f_conv are inferred from Fig. 1 of Oh & Kroupa (2016).

2.6 Age uncertainty

In this work we assumed that stars within a cluster have the same age. Observationally inferred cluster ages, however, can carry uncertainties, which propagate into the inferred stellar evolutionary state and (in our framework) into both the expected dynamical-ejection fraction and the Hα luminosity. J23 discuss age uncertainties at the level of ≈1–2 Myr.

To account for this, we treated the age uncertainty as a sensitivity test by re-evaluating the stellar masses implied by a given ionising output at a shifted age. Concretely, we kept the previously computed L_Hα for each star and inferred the corresponding initial mass (m_ini) using the L_Hα − m_ini relation at the shifted age. We then used the evolutionary tracks (Sect. 2.3) to obtain the stellar mass at that age. When the mapping is not single-valued at older ages (multiple m_ini yielding the same L_Hα), we chose the solution closest to the original mass to minimise discontinuous jumps. If L_Hα exceeded the maximum value produced at the shifted age in our precomputed grid, we treated this as an indication that the star must be younger than the shifted age within our modelling framework and therefore used the youngest age bin considered for the L_Hα − m_ini mapping.

2.7 Contamination by field stars

The pixel scale of the camera used by J 23 is 0″.0396 pixel⁻¹. J23 chose an aperture radius of 5 pixels, corresponding to 8–10 pc for the distances (4.3–5.1 Mpc) of the selected galaxies (Jung et al. 2023). The targeted clusters have diameters smaller than 9 pc, and so crowded environments and aperture-based selection can lead to residual local contamination and imperfect membership assignment.

We accounted for this using a controlled sensitivity test and adding additional stars to each model cluster. The added stars have masses m ∈ (0.1, 150) M_⊙ and were drawn from the cluster’s IMF up to a total added stellar mass of 10% of the cluster mass. Each added star was assigned a random age between 1 and 6 Myr. We computed the corresponding L_Hα as described in Sect. 2.4 and added it to the cluster total, together with the corresponding mass contribution. We adopted the 10% contamination level and the 1–6 Myr age range as a conservative upper-limit bracketing case, intended to test how strongly residual contamination could affect the inferred cluster distributions.

2.8 Mean and standard deviation

To quantify the results better, the weighted mean is calculated with Eq. (6) as J23 did. The standard deviation of the mean is determined with Eq. (7) from Taylor (1997), where N describes the number of star clusters: $$⟨\frac{L_{\mathrm{H}\alpha }}{M_{\mathrm{c}\mathrm{l}}}⟩={\log }_{10}\left(\frac{\sum _i{L}_{\mathrm{H}\alpha ,i}}{\sum _i{M}_{\mathrm{c}\mathrm{l},i}}\right)$$(6) $${\sigma }_{⟨L_{\mathrm{H}\alpha }/M_{\mathrm{c}\mathrm{l}}⟩}=\frac{1}{\sqrt{N}}\sqrt{\frac{1}{N}\sum _i{({\log }_{10}\left(\frac{L_{\mathrm{H}\alpha ,i}}{M_{\mathrm{c}\mathrm{l},i}}\right)-\frac{1}{N}\sum _i{\log }_{10}\left(\frac{L_{\mathrm{H}\alpha ,i}}{M_{\mathrm{c}\mathrm{l},i}}\right))}^2}.$$(7)

J23 divided the clusters into three different mass ranges, namely (5.17–23.17) × 10² M_⊙, (2.35–14.52) × 10³ M_⊙ and (1.881–8.567) × 10⁴ M_⊙. The clusters of this work were divided into similar bins but not exactly the same ones, namely (5.0–22.9) × 10² M_⊙, (2.3–14.9) × 10³ M_⊙, and (1.5–10.0) × 10⁴ M_⊙ because of different cluster masses. The three mass ranges are called Bins 1, 2, and 3, respectively. In addition, the weighted mean and standard deviation were calculated for all clusters (Bin 5) of this work, so for masses in (316–10⁵) M_⊙ and for clusters that are in the range of the clusters observed by J23 (Bin 4). For a better comparison to J23, the same number of clusters as in J23 is picked randomly from the simulated data and the mean is calculated for the respective bin. This is repeated 1000 times from which the median and the median absolute deviation (MAD) are determined. Note that throughout this work, we drop the units for conciseness, but solar units are always implied: $⟨L_{\mathrm{H}\alpha }/M_{\mathrm{c}\mathrm{l}}⟩=⟨\frac{L_{\mathrm{H}\alpha }/M_{\mathrm{c}\mathrm{l}}}{L_{\mathrm{H}\alpha ,\odot }/M_{\odot }}⟩$.

3 Results

In Sects. 3.1 and 3.2 the computational results of the stellar evolution and the Hα luminosity are presented. In Sect. 3.3 the sampled clusters are depicted under different assumptions, such as dynamical ejections and observational uncertainties.

3.1 Stellar evolution

In Fig. 1 the stellar mass evolution of stars with ages of 1–4 Myr and different metallicities is displayed. The stellar mass evolution for stars with ages of 5–7 Myr is displayed in Fig. C.1.

For low m_ini, the mass at an age of 1 Myr does not change but the mass of stars with m_ini > 60 M_⊙ decreases. The higher the initial mass, the stronger the mass loss. The stellar mass evolution is only displayed for m_ini ≤ 150 M_⊙ because this is the upper mass limit of the IMF and no mergers leading to higher masses are included in this analysis. At an age of 2 Myr, stars with m_ini > 40 M_⊙ undergo notable mass loss, especially stars with solar metallicity. At an age of 3 Myr, no stars heavier than m_ini ≈ 140 M_⊙ exist. Furthermore, the dependence between m_ini and m becomes more complex. At an age of 4 Myr, stars up to a mass of 30 M_⊙ have not changed their mass appreciably, while heavier stars can lose more than half of their initial mass. Stars with m_ini > 65 M_⊙ have lost their entire mass and are therefore no longer present at this age. The higher the metallicity, the stronger the mass loss. For low-mass stars, the metallicity does not change the evolution significantly at the depicted ages. Due to mass loss, a star with a higher initial mass can become less massive than a star with a lower initial mass at a given age. For example, a star with m_ini ≈ 140 M_⊙ can reach a similar mass at 3 Myr as a star with m_ini ≈ 30 M_⊙.

Fig. 1 Stellar mass evolution for ages of 1–4 Myr obtained as described in Sect. 2.3. The mass of a star at the respective age, m, is plotted against the initial mass, mini. The colour indicates the metallicity, and the marker the age. The black line displays the case m = mini.

3.2 Ho luminosity

In Fig. 3 the dependence of the Hα luminosity, L_Hα, on the initial stellar mass, m_ini, is displayed for different metallicities and ages of 1–4 Myr. These results are obtained with Pégase as explained in Sect. 2.4. The plots of the Hα luminosity against the initial stellar mass split by age (1–7 Myr) are shown in Fig. A.1.

While in the higher-mass regime (m > 25 M_⊙), the curve is flat, for lower initial masses, L_Hα increases strongly. Thus, a different mass up to m_ini ≈ 25 M_⊙ has a higher impact on L_Hα than a change of mass in the higher-mass regime. The data points for stars with an age of 3 Myr are mostly covered by the data points of stars with an age of 4 Myr. At these ages, for m_ini > 70 M_⊙, the values change without apparent correlation. Since the model’s upper mass limit is 120 M_⊙, L_Hα is extrapolated up to m_max = 150 M_⊙ for ages at which these stars are still alive. There is no visible correlation between L_Hα and metallicity.

3.3 Sampled clusters

With the sampling methods explained in Sect. 2.2, 160 star clusters with M_ecl = 10^2.5−10⁵ M_⊙ and metallicities of Z = {0.014, 0.004, 0.008} are sampled (see Sect. 2 for details on the model cluster-mass selection and the matching procedure used for a fair comparison to J23). The resulting clusters, for which only the stellar evolution is accounted for, are displayed in Fig. 2 and described in Sect. 3.3.1. The clusters on which dynamical ejections are applied are shown in Sect. 3.3.2. The results in which the age uncertainty and the field stars are included are depicted in Sects. 3.3.3 and 3.3.4, respectively.

In the upper panels of Figs. 2–6, the mass of the most massive star, m_max, is plotted against the mass of the stars in the cluster, M_cl at different ages. On the left, the optimally sampled clusters are displayed and on the right, the randomly sampled clusters. In the lower panels, the ratio of the Hα luminosity and M_cl, L_Hα/M_cl, in dependence of M_cl is depicted.

Fig. 2 Optimally and randomly sampled clusters with ages of 1–4 Myr. The colour indicates the metallicity of the cluster, and the marker the age. Top left: mass of the most massive star in the cluster, mmax, against the stellar mass of the cluster, Mcl, for optimally sampled clusters. Top right: same as the top left but for randomly sampled clusters. Bottom left: ratio of the cluster’s Hα luminosity and Mcl against Mcl for optimally sampled clusters. Bottom right: same as the bottom left but for randomly sampled clusters.

Fig. 3 Hα luminosity (logarithmic scale) of stars with ages of 1–4 Myr as a function of the initial mass, mini. The colour indicates the metallicity, and the marker indicates the age. Note that the 1–3 Myr points are largely superposed by the 4 Myr data. In Appendix A, the data for each age are displayed in a separate figure for clarity.

3.3.1 Including stellar evolution

For the optimally sampled clusters displayed in the upper left panel of Fig. 2, a clear correlation between m_max and M_cl, the m_max–M_ecl relation, is visible. Different plotted ages lead to changing break points at which the relation flattens. During a cluster’s evolution, the stars develop and experience mass loss, which increases with initial stellar mass (see Fig. 1). So, when the cluster ages, the mass of the most massive star decreases, which results in the lower break points. For the oldest clusters, some deviations from the m_max–M_ecl relation are visible. This is caused by the further evolution of the stars as displayed in Fig. 1.

A lower metallicity results in a more massive m_max. This is a result of the varying IMF becoming top-heavier with lower metallicities (see Sect. 2.1).

The randomly sampled clusters do not show a clear relation between m_max and M_cl. However, the samples exhibit plateaus that vary with age and metallicity due to the upper-mass limit m_up = 150 M_⊙. Some clusters are missing or do not follow the respective plateau, especially clusters with M_c < 1000 M_⊙. Fewer of these clusters host a star with a mass close to the upper mass limit because of the lower total mass. Furthermore, the randomly sampled clusters are able to host heavier stars and a lower number of objects than the optimally sampled clusters. After the death of the heaviest star, the second most massive star is the new most massive star whose mass is different to the original m_max, which leads to a plateau. An important aspect is that clusters with an age of 4 Myr are shifted towards lower M_cl. In the random-sampling case, some clusters can consist of only a few very massive stars; once these stars evolve off, the cluster loses a substantial fraction of its stellar mass and shifts to lower M_cl.

In the lower panels, we plot the ratio L_Hα/M_cl as a function of M_cl. For the optimally sampled clusters, L_Hα/M_cl increases with increasing M_cl, indicating that the ionising output grows faster than the stellar mass. This is expected because clusters with M_cl ≳ 500 M_⊙ typically contain stars more massive than ≈20 M_⊙, for which L_Hα is orders of magnitude higher than for lower-mass stars. The slope of the L_Hα/M_cl–M_cl relation decreases towards higher cluster masses, consistent with Fig. 3, because the Hα luminosity of individual stars changes only weakly for m_ini ≳ 25 M_⊙ in our models. Even once m_max approaches the upper-mass limit, L_Hα/M_cl can still increase because additional ionising sources are added, rather than a single more luminous star. For clusters older than ≈3 Myr, stellar mass loss and the more complex mapping between L_Hα and m_ini weaken the correlation. Finally, we do not find a strong systematic metallicity trend in L_Hα/M_cl over the range considered here.

The ratio L_Hα/M_cl of the randomly sampled clusters does not show an apparent correlation to M_cl. Here, no gradient is visible compared to the optimally sampled clusters.

In J23, most clusters lie in the first bin (see Table D.1). In our model set, the optimally and randomly sampled clusters are distributed more uniformly across the first three bins. The fourth bin indicates that not all sampled clusters fall within the mass range analysed by J23, because our initial model grid spans a wider range to follow the evolution with time. The medians of the weighted means of the randomly sampled clusters have a range of ⟨L_Hα/M_cl⟩ = (35.55–35.65) in the first three bins and the optimally sampled clusters of ⟨L_Hα/M_cl⟩ = (34.31–35.23). The values of the weighted mean in J23 range from ⟨L_Hα/M_cl⟩ = (33.900–34.176), which are lower than the results from both sampling methods.

Fig. 4 Same as Fig. 2 but additionally including dynamical ejections with a 2σ higher ejection fraction.

3.3.2 Including dynamical ejections

The results for an ejection fraction within 2σ confidence are plotted in Fig. 4. With the higher ejection fraction, the m_max–M_ecl relation for the optimally sampled clusters becomes more dispersed. For clusters with M_cl ≲ 1500 M_⊙ it is no longer clearly visible as a tight relation, but remains as an upper envelope. The ejection of the most massive stars reduces m_max, while M_cl changes less because most of the cluster mass resides in low-mass stars that are unlikely to be ejected. As a result, the clusters move downwards and slightly leftwards in the diagram, but do not show the extreme shifts towards very low M_cl seen in some randomly sampled cases. The high-mass clusters remain largely unchanged because the adopted ejection-fraction model is limited to (and decreases towards) the highest masses covered by the simulations. In the lower panel, the ratios of the optimally sampled clusters appear more scattered in the range where dynamical ejections are applied. Older clusters tend to show lower L_Hα/M_cl than younger clusters. Since M_cl decreases, L_Hα must decrease even more strongly, consistent with the preferential ejection of the most massive (ionising) stars.

One cluster lies well below the main population (left panels at log₁₀(M_cl/M_⊙) ≈ 2.5). The cluster has lost its most massive star but has not shifted markedly to lower M_cl, implying that most of its remaining mass resides in low-mass stars. Since these stars contribute little to the ionising photon budget, the loss of a massive star reduces L_Hα substantially and therefore decreases L_Hα/M_cl.

In the upper-right panel of Fig. 4, where the randomly sampled clusters are displayed, a gap appears at log₁₀(M_cl/M_⊙) ≈ 3.0 after applying ejections. This can occur because some randomly sampled clusters contain a small number of massive stars and comparatively few low-mass stars; when the massive stars are ejected, the remaining stellar mass drops sharply and the cluster shifts leftwards in the diagram. This effect leads to a correlated looking population of clusters for M_cl ≈ (10–100) M_⊙. Note that the clusters populating the lower-left corner consist of only one or a few remaining stars. Furthermore, the randomly sampled clusters are overall shifted further to the left because of the strong mass loss of the most massive stars and the loss of stars due to dynamical ejections.

In the lower-right panel, the stars are not as dispersed as in the previous plots and are mostly populated in the range log₁₀(L_Hα/M_cl) ≈ (35.0–36.5) except for the clusters that have lost a significant amount of their mass, which are populated lower and on the left. Furthermore, clusters with masses of around 100 M_⊙ have higher values of L_Hα/M_cl now. This indicates that M_cl has decreased significantly with some stars, which led to a high L_Hα still remaining in the cluster.

With the higher ejection fraction, the randomly sampled clusters have lost a significant amount of mass and 20 clusters moved out of the mass range by J23 to lower masses (see Table D.3). The values of ⟨L_Hα/M_cl⟩ are a little bit lower than before but still not in agreement with J23.

Fig. 5 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 1 Myr.

Fig. 6 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction, an age uncertainty of 1 Myr, and contamination of field stars.

3.3.3 Including an age uncertainty

J23 question the age of the clusters and estimate that the clusters might be 1–2 Myr older. In our Figs. 5, D.3, and D.4, the clusters are assumed to be 1 Myr, 2 Myr, and 3 Myr older, respectively.

When we include an age uncertainty, the dispersion of the optimally sampled clusters increases in the upper-left panels. As the age uncertainty increases, the plateaus merge more strongly for both sampling methods. Older clusters tend to shift upwards because stellar masses are reassigned based on L_Hα. The change in stellar mass increases with age (Figs. 1 and C.1), and the maximum L_Hα produced by the L_Hα–m_ini mapping decreases at older ages. As described in Sect. 2.6, if L_Hα exceeds the maximum value available at the shifted age in our precomputed grid, we used the youngest age bin considered for the L_Hα– m_ini mapping.

In the lower panel of Fig. 5, the clusters also appear randomly populated except for the higher-mass clusters because of the lower dynamical ejection fractions for massive clusters. For larger assumed age uncertainties, the older star clusters move upwards, indicating that the inferred cluster mass decreases with increasing age uncertainty in this prescription.

The randomly sampled clusters show plateaus as before, but some clusters at ages of 3–4 Myr populate the upper plateau. This behaviour arises from the age-uncertainty prescription, which reassigns stellar masses based on a shifted-age L_Hα–m_ini mapping (Sect. 2.6). Furthermore, about 25 low-mass clusters develop an apparent m_max–M_cl correlation under these assumptions.

In the lower-right panel, the clusters are mostly populated at values between log₁₀(L_Hα/M_cl) ≈ 35–36 except some older clusters that build a similar correlation like in the upper panel. Also, the spread of the randomly sampled clusters in the lower panel does not change significantly with age uncertainty.

Compared to before, the number of randomly sampled clusters in the second and third bin decreased. On the other hand, the total number of randomly sampled clusters has not changed much, which implies that the clusters moved out of the mass range set in J23. Furthermore, the weighted mean of the randomly sampled clusters increased while the one of the optimally sampled clusters decreased as displayed in Table D.4. Since higher age uncertainties do not change the lower panels except for the massive optimally sampled clusters, which are not included in the set of clusters in J23, an age uncertainty of 1 Myr is assumed in the following plots.

Fig. 7 Left: sum of the lower-left panels of Figs. 5 and 6. All ages are represented by the dots. The colour indicates the metallicity as before. The black box displays the limits from J23, and the red grid shows the bins used for the χ2 test (Sect. 4.3). Right: same as the left but for randomly sampled clusters. The grey data points represent the observational data from J23 (see their Fig. 4 for error bars, which are omitted here for clarity).

3.3.4 Including field stars

As described in Sect. 2.7, stars may be mis-assigned to clusters. We accounted for this uncertainty by adding contaminating stars drawn from the same IMF as the cluster (Sect. 2.7).

With added field-star contamination, the plateaus become less distinct, as shown in the upper right panel of Fig. 6. Furthermore, all clusters shift upwards. Because additional stars are added, some clusters can contain a new m_max. In addition, M_cl increases, which shifts all clusters slightly to the right.

In the lower panel, clusters with M_cl ≈ 1000 M_⊙ are more dispersed. Some are shifted significantly upwards because of a new most massive star that leads to a higher L_Hα whose increase must exceed the one of M_cl because of the higher ratio L_Hα/M_cl. High-mass clusters remain populated mostly in the same region, consistent with the upper panel. This is because high-mass clusters already host a larger number of massive stars than low-mass clusters, so adding a small contaminating component typically does not change their configuration as strongly.

The randomly sampled clusters change less under this prescription. Since they are already stochastically drawn, adding additional IMF-drawn stars primarily increases M_cl and produces a modest shift towards higher M_cl.

Referring to Table D.5, adding field stars does not substantially change the resulting number of star clusters, weighted mean, and standard deviation for the randomly sampled clusters. For the optimally sampled clusters, ⟨L_Hα/M_cl⟩ increases compared to the case without field stars, consistent with Fig. 6.

4 Discussion

In this section the results from Sect. 3 are discussed. First the influence of field stars and dynamical ejections are investigated. In Sect. 4.2 the final plot (Fig. 7) and values are compared to the results by J23 and a χ² test is applied in Sect. 4.3.

4.1 Dynamical ejections and field stars

In Fig. D.2 the clusters are displayed with dynamical ejections with a mean ejection fraction, an age uncertainty of 1 Myr, and randomly added field stars. Compared to Fig. 6, clusters with M_cl ≈ 1000 M_⊙ are spread along the y-axis. On the lower right (log₁₀(M_cl/M_⊙) > 4), this distribution is missing caused by the field stars. Thus, dynamical ejections move the star clusters downwards because of the ejections of massive stars, which reduce L_Hα more than M_cl. Adding field stars moves the star clusters upwards in the plot because of the higher ratio of massive stars in the star cluster than before. This affects optimally sampled star clusters with masses around 1000 M_⊙ the most. Furthermore, the dynamical ejections lead to the population of randomly sampled clusters in the lower-left corner, which is reduced with the added field stars.

The criterion for a star to be ejected is its distance from the cluster centre to be greater than 3 × r_h at 3 Myr, and its velocity to be larger than the escape velocity at its distance (Oh et al. 2015). Star clusters expand in their evolution up to ten times of their initial half-mass radius (Banerjee & Kroupa 2017). Thus, the clusters of this paper would have r_h ≈ 3 pc at an age of 3 Myr. This corresponds to a minimum distance of 9 pc. Furthermore, the ejected systems have generally high velocities, which make them to move away from the cluster fast (Oh et al. 2015). Thus, those stars would no longer be in the frame of the aperture. A last note is that 3 × r_h is only the minimum distance; the real distance of the ejected systems might be greater.

4.2 Final comparison

Because the level of residual contamination and mis-assigned stars is uncertain, we treated the ’field star’ prescription as a bracketing case. In Fig. 7 we therefore combine the model realisations with and without added field stars to illustrate the envelope of plausible outcomes once dynamical ejections are included. The corresponding summary statistics are listed in Table D.6. For the comparison to J23, we restricted the model to the mass range covered by the observed sample and down-sampled the highest-mass model clusters to reflect the scarcity of high-mass clusters in J23.

To better quantify the results, the weighted mean and standard deviation are calculated as described in Sect. 2.8. The resulting values for clusters on which dynamical ejections with a 2σ higher ejection fraction, an age uncertainty of 1 Myr and on one set field stars and on the other no field stars are applied are listed in Table D.6. The resulting weighted means of the optimally sampled clusters are in better agreement with J23 for all bins. The randomly sampled clusters have too high means. This is also reflected in Fig. 7 where most of the randomly sampled clusters lie outside the range of J23.

The standard deviation provides a measure of dispersion in the data points. Except for the last bin, σ_{⟨LHα/Mcl⟩} for the optimally sampled clusters is about twice as high as for the randomly sampled clusters, indicating a larger dispersion. This is consistent with Fig. 7, where the randomly sampled clusters lie mainly in log₁₀(L_Hα/M_cl) ≈ 35–36 for M_cl > 100 M_⊙, while the optimally sampled clusters span approximately log₁₀(L_Hα/M_cl) ≈ 32–35. In Bin 3, σ_{⟨LHα/Mcl⟩} of the randomly sampled clusters is higher than for the optimally sampled clusters. That can be explained with the missing dynamical ejections in this range. However, these clusters do not influence the final comparison in Fig. 7 because J23 do not observe a large number of massive star clusters. Therefore, the number of clusters is too low for a reliable comparison. Furthermore, σ_{⟨LHα/Mcl⟩} cannot be compared in a strictly like-for-like way to the values reported by J23 because the exact estimator used there is not fully specified. If their dispersion measure is sensitive to sample size, then differences in the number of clusters per bin can also affect the reported values. We therefore interpret the dispersion comparison with caution and focus on the overall distributional agreement in the diagnostic plots.

Referring to Table D.6, the number of optimally sampled clusters does not change in the respective mass bins significantly, which indicates that M_cl does not change with the added influences. On the other hand, the number of randomly sampled clusters decreases and M_cl shrinks because of the mass loss caused by the added influences. Furthermore, the added effects lead to a higher σ_{⟨LHα/Mcl⟩} for the first four bins of the optimally sampled clusters while the values of the randomly sampled clusters decrease.

To facilitate a direct visual comparison to J23, we indicate the observed range from J23 with a black box in Fig. 7. We down-sampled the highest-mass model clusters so that the comparison is not dominated by cluster masses that are sparsely represented in the observed sample.

Referring to Fig. 7, the randomly sampled clusters are not in the range of the clusters analysed by J23, except for a few. Furthermore, in the cluster mass range, the ratios of L_Hα and M_cl do not vary that much as in J23. For the optimally sampled clusters, displayed in the left panel, the clusters are spread similar to J23 and are in the same range as their results. So, the optimally sampled clusters are more dispersed than the randomly sampled ones. This may appear counter-intuitive at first glance, but it illustrates that an initially ordered (deterministic) stellar population can evolve into a distribution that looks ’stochastic’ once evolutionary processes and observational uncertainties are taken into account.

4.3 χ² test

To compare the distributions further, we applied a χ² test. The bin edges were chosen to be ${\log }_{10}\left(\frac{L_{\mathrm{H}\alpha }/M_{\mathrm{c}\mathrm{l}}}{L_{\odot }/M_{\odot }}\right)=32,33,34,35$ and log₁₀(M_cl/M_⊙) = 2.5, 3.0, 3.5, 4.0 and are displayed in Fig. 7. The number of data points for each bin cell was determined for the distribution in J23, n_J,i, and for the calculated distribution, n_c,i. The total number of data points is $N_{\mathrm{J}}={\sum }_i^N{n}_{\mathrm{J},\mathrm{i}}$ for the data in J23 and $N_{\mathrm{c}}={\sum }_i^N{n}_{\mathrm{c},\mathrm{i}}$ for the calculated data. In accordance with Eq. (14.3.3) in Press et al. (1986), the χ² of binned data with an unequal number of data points was calculated as ${\chi }^2={\sum }_i^N\frac{{(\sqrt{N_{\mathrm{J}}/N_{\mathrm{c}}}{n}_{\mathrm{c},\mathrm{i}}-\sqrt{N_{\mathrm{c}}/N_{\mathrm{J}}}{n}_{\mathrm{J},\mathrm{i}})}^2}{n_{\mathrm{c},\mathrm{i}}+n_{\mathrm{J},\mathrm{i}}}$ and divided by the degrees of freedom, d.o.f = N_bin − 1, where N_bin is the number of bins. For the optimally sampled clusters, $\frac{{\chi }^2}{\text{d.o.f}}\approx 4$, and for the randomly sampled clusters, $\frac{{\chi }^2}{\text{d.o.f}}\approx 17$. The corresponding p-values are ≈10⁻⁴ and ≈10⁻²⁵, which correspond to the 4 and 11σ regions. This underlines the disagreement between the observed and randomly sampled clusters.

4.4 Further comments

There are additional effects, such as gas expulsion (Kroupa et al. 2001; Dinnbier & Kroupa 2020; Wirth et al. 2024), which we did not include in this analysis. Gas expulsion can reduce the bound stellar mass and potentially introduce additional dispersion in the observable relations (Weidner et al. 2007). However, the optimally sampled models already show substantially better agreement with the observational distribution than the purely randomly sampled models under the set of effects explored here. In particular, the combination of stellar evolution, dynamical ejections, age uncertainty, and contamination can broaden the distribution of optimally sampled clusters in Fig. 7, whereas the randomly sampled clusters remain concentrated in a narrower region. We therefore conclude that evolutionary and observational effects can make an initially regulated stellar population appear more stochastic, and that a purely random sampling scenario does not reproduce the observed distribution within the framework investigated here.

5 Conclusion

In this study we analysed the evolution of star clusters using both optimal and random sampling methods. Clusters were modelled with masses between 10^2.5 and 10⁵ M_⊙ and metallicities of Z = {0.014, 0.004, 0.008}. We examined how these sampling methods influence the m_max–M_ecl relation and the ratio L_Hα/M_cl.

This study shows that, within the modelling framework explored here, optimal sampling yields distributions that are consistent with the extragalactic star cluster observations considered, while purely random sampling produces no agreement. Our modelling further suggests that the observed scatter in the LEGUS clusters can be accounted for by a combination of stellar evolution, dynamical ejections, age uncertainties, and residual contamination by field stars, rather than by the absence of the m_max–M_ecl relation.

For reproducibility, we provide our Python-based random-sampling implementation (linked in Sect. 2.2.1). Overall, our results support a highly self-regulated interpretation of star formation and are consistent with key aspects of the integrated galaxy-wide IMF framework.

Acknowledgements

We acknowledge support through the DAAD-EasternEurope Exchange grant at Bonn University and corresponding support from Charles University. P.K. acknowledges support through grant 26-217745 from the Czech Grant Agency. T.J. acknowledges the support from the MUNI Award in Science and Humanities (MUNI/I/1762/2023).

References

Appendix A Hα luminosity

The downloaded files from https://www2.iap.fr/users/fioc/PEGASE.html comprise the stellar evolution tracks of the following metallicities Z = {0.0001, 0.0004, 0.004, 0.008, 0.02, 0.05, 0.1}. Furthermore, various IMFs are available, where the Kroupa et al. (1993) IMF is chosen as IMF input parameter. Kroupa et al. (1993) describe the IMF as a three-part power law, where α₃ = 2.7 for stars more massive than 1 M_⊙, α₂ = 2.2 in the mass range 0.5 ≤ m/M_⊙ ≤ 1 and α₁ = 1.3 in the range 0.1 < m/M_⊙ ≤ 0.5. In addition, multiple filter and calibration options are available, which are not important for this work.

The first step in the execution is to create SSPs. As boundary masses, the masses 1 M_⊙ below and above the concerned mass are chosen. For example, if a SSP for m_star = 10 M_⊙ is set up, the lower boundary mass is 9 M_⊙ and the upper 11 M_⊙. For the supernovae ejecta, model B is selected, which contains the calculations of Woosley & Weaver (1995). Furthermore, ejecta due to stellar winds in high-mass stars are accounted for.

After generating the SSPs, the scenarios are created. The fraction of close binary systems is set to 0 to allow normalisation to a single star later. The initial metallicity of the ISM is set to the respective metallicity of the galaxy or to the solar metallicity. Furthermore, the galaxy is set to be already constituted and the star formation scenario to be an instantaneous burst SFR(t) = δ(t). The fraction (in mass) of the star formation rate used to form substellar objects is set to 0. In addition, the galactic winds and the extinction are set to 0, while the nebular emission is accounted for. After executing spectra, the normalised luminosities of the line with a wavelength of λ = 6563 Å, which is similar to the wavelength of the Hα line, are extracted.

Fig. A.1 Same as Fig. 3 but split by age for 1–7 Myr (in order from 1 Myr to 7 Myr).

Appendix B Dynamical ejections

Appendix C Stellar evolution

Fig. C.1 Same as Fig. 1 but for ages of 5 Myr (left), 6 Myr (middle), and 7 Myr (right).

Appendix D Sampled clusters

In Fig. D.1 the results after applying dynamical ejections with the mean ejection fraction are plotted.

In comparison to Fig. 2, the m_max–M_ecl relation in the upper-left panel of Fig. D.1 appears more dispersed but remains visible. In particular, clusters older than 2 Myr with M_cl ≈ 900–1700 M_⊙ show the largest changes relative to the initial plot. This is consistent with the peak of the ejection fraction at this cluster-mass scale (see Fig. 2 in Oh et al. 2015). A higher ejection probability for the most massive star reduces m_max and, to a lesser extent, M_cl, shifting clusters downwards and slightly leftwards. Note that the lower panels use different y-axis limits than in Fig. 2, which affects the visual impression. For clusters with M_cl ≳ 10^4.5 M_⊙, we did not apply dynamical ejections because the fraction of ejected massive stars becomes negligible (Oh et al. 2015). We therefore leave the highest-mass model clusters unchanged with respect to dynamical ejections.

The randomly sampled clusters change less, except that more objects populate the low-M_cl regime. These clusters are typically at 4 Myr and have undergone substantial mass loss and dynamical ejections, which reduces the remaining stellar mass reservoir. Optimally sampled clusters are less affected by stellar evolution and dynamical ejections because their stellar populations contain a larger fraction of low-mass stars, which are less likely to be ejected and lose less mass over these young ages.

The resulting values of the number of clusters, weighted mean and standard deviation are listed in Table D.2. Compared to before, the mass of the randomly sampled clusters is reduced because there are more clusters in the first bin than before and less in the second and third bins. Furthermore, one cluster is dissolved or has M_cl < 10 M_⊙ because there are only 159 clusters in the fifth bin. The distribution of optimally sampled clusters does not change significantly. Only one cluster moved from the third bin into the second. This supports the argumentation from before that the M_cl of the randomly sampled clusters is more sensitive to dynamical ejections because of the higher amount of massive stars. The values of ⟨L_Hα/M_cl⟩ and σ_{⟨LHα/Mcl⟩} are still in the same range as before for the randomly and optimally sampled clusters.

Fig. D.1 Same as Fig. 2 but including dynamical ejections with the mean ejection fraction.

Fig. D.2 Same as Fig. 2 but including dynamical ejections with a mean ejection fraction, an age uncertainty of 1 Myr, and field stars.

Fig. D.3 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 2 Myr.

Fig. D.4 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 3 Myr.

All Tables

All Figures

	Fig. 1 Stellar mass evolution for ages of 1–4 Myr obtained as described in Sect. 2.3. The mass of a star at the respective age, m, is plotted against the initial mass, mini. The colour indicates the metallicity, and the marker the age. The black line displays the case m = mini.
In the text

Fig. 2 Optimally and randomly sampled clusters with ages of 1–4 Myr. The colour indicates the metallicity of the cluster, and the marker the age. Top left: mass of the most massive star in the cluster, mmax, against the stellar mass of the cluster, Mcl, for optimally sampled clusters. Top right: same as the top left but for randomly sampled clusters. Bottom left: ratio of the cluster’s Hα luminosity and Mcl against Mcl for optimally sampled clusters. Bottom right: same as the bottom left but for randomly sampled clusters.

In the text

	Fig. 3 Hα luminosity (logarithmic scale) of stars with ages of 1–4 Myr as a function of the initial mass, mini. The colour indicates the metallicity, and the marker indicates the age. Note that the 1–3 Myr points are largely superposed by the 4 Myr data. In Appendix A, the data for each age are displayed in a separate figure for clarity.
In the text

	Fig. 4 Same as Fig. 2 but additionally including dynamical ejections with a 2σ higher ejection fraction.
In the text

	Fig. 5 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 1 Myr.
In the text

	Fig. 6 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction, an age uncertainty of 1 Myr, and contamination of field stars.
In the text

Fig. 7 Left: sum of the lower-left panels of Figs. 5 and 6. All ages are represented by the dots. The colour indicates the metallicity as before. The black box displays the limits from J23, and the red grid shows the bins used for the χ2 test (Sect. 4.3). Right: same as the left but for randomly sampled clusters. The grey data points represent the observational data from J23 (see their Fig. 4 for error bars, which are omitted here for clarity).

In the text

	Fig. A.1 Same as Fig. 3 but split by age for 1–7 Myr (in order from 1 Myr to 7 Myr).
In the text

	Fig. C.1 Same as Fig. 1 but for ages of 5 Myr (left), 6 Myr (middle), and 7 Myr (right).
In the text

	Fig. D.1 Same as Fig. 2 but including dynamical ejections with the mean ejection fraction.
In the text

	Fig. D.2 Same as Fig. 2 but including dynamical ejections with a mean ejection fraction, an age uncertainty of 1 Myr, and field stars.
In the text

	Fig. D.3 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 2 Myr.
In the text

	Fig. D.4 Same as Fig. 2 but including dynamical ejections with a 2σ higher ejection fraction and an age uncertainty of 3 Myr.
In the text