2.1 Praesepe

For Praesepe, the centre was adopted as $\alpha=08^{\text{h}}40^{\text{m}}13^{\text{s}}$ , $\delta=19^{\circ}37'16''$ (Tarricq et al. Reference Tarricq2021). Only Gaia DR3 sources within $15^{\circ}$ of this position were considered in the analysis, which corresponds to radii of $\approx$ 36–50 pc at the distances of the three studied clusters. These radii are much larger than the typical tidal radius of open clusters of $\approx6$ pc (Piskunov et al. Reference Piskunov, Schilbach, Kharchenko, Röser and Scholz2008) and thus should include all bound cluster members. Member candidates were then distinguished from background stars using their proper motions, parallaxes, and photometry.

A $\chi^2$ test was used to filter member stars based on proper motions and parallaxes. More specifically, a given star i has a $\chi^2$ value given by

(1) \begin{equation} \chi_i^2= \frac{\left(\omega_i-\bar{\omega}\right)^2}{e_{\omega,i}^2 + \sigma_\omega^2} + \frac{\left(\unicode{x03BC}_{\alpha,i} - \bar{\unicode{x03BC}}_\alpha\right)^2}{e_{\alpha,i}^2 + \sigma_\alpha^2} + \frac{\left(\unicode{x03BC}_{\delta,i} - \bar{\unicode{x03BC}}_\delta\right)^2}{e_{\delta,i}^2 + \sigma_\delta^2} \end{equation}

and passes the initial filtering if

(2) \begin{equation} \chi^2_i \leq \chi^2_{0.99,\nu=3}\approx11.345\end{equation}

where $\omega_i$ and $\bar{\omega}$ are the parallax of star i and the mean cluster parallax, $\unicode{x03BC}_{\alpha,i}$ and $\bar{\unicode{x03BC}}_\alpha$ are the proper motion of star i and the mean cluster proper motion in right ascension, $\unicode{x03BC}_{\delta,i}$ and $\bar{\unicode{x03BC}}_\delta$ are the proper motion of star i and the mean cluster proper motion in declination, $e_{\omega,i}$ , $e_{\alpha,i}$ and $e_{\delta,i}$ are the errors of each measurement, $\sigma_\omega$ is the intrinsic spread of parallaxes due to the line of sight depth of the cluster, $\sigma_\alpha$ and $\sigma_\delta$ are the internal 1D velocity dispersions and $\nu$ is the number of degrees of freedom.

For the internal velocity dispersions, we make conservative estimates of $\sigma_\alpha=\sigma_\delta=1\,\mathrm{mas\,yr^{-1}}$ which correspond to velocities of $\approx0.9\,\mathrm{km\,s^{-1}}$ at the distance of Praesepe. This value is similar to typical estimates of the intrinsic velocity dispersion in open clusters, which are generally of order $\lesssim1\,\mathrm{km\,s^{-1}}$ , and was chosen to ensure that the selection encompasses both the intrinsic kinematic spread of cluster members and the contribution from measurement uncertainties in Gaia proper motions. The parallax dispersion was chosen to be $\sigma_\omega=0.25$ , which corresponds to a line of sight radius of $\approx8$ pc at the distance of Praesepe, comparable to or slightly larger than the expected physical extent of the cluster.

The membership threshold is also a conservative value, using the 99th percentile with 3 degrees of freedom of a $\chi^2$ distribution, meaning that only 1% of true member stars are likely to be excluded due to measurement error assuming Gaussian errors. While these values will ensure nearly all member stars are detected, the chance of including contaminating background stars in the data set is also increased. Any remaining contamination can subsequently be filtered out using photometry measurements and cross referencing with additional catalogues.

For the cluster mean parallax and proper motion, we began with the values $\bar{\omega}=5.371$ mas, $\bar{\unicode{x03BC}}_\alpha=-36.047\,\mathrm{mas\,yr^{-1}}$ and $\bar{\unicode{x03BC}}_\delta=-12.917\,\mathrm{mas\,yr^{-1}}$ from Gaia DR2 (Gaia Collaboration et al. Reference Collaboration2018). We then implemented a routine where member stars are selected using these initial reference values using the $\chi^2$ filtering. The new mean values for these three values, calculated from the resulting member population, are then used in the next iteration of the filtering on the whole population again. This process was repeated until convergence is reached and the selected member stars are consistent between successive iterations. The result of this routine can be seen in proper motion space in Figure 1, where we detect a compact over-density in proper motion and parallax out from the background stellar population.

Figure 1.

Proper motion distribution of stars in Praesepe in right ascension and declination. Black points show background stars whereas red points indicate cluster members identified through PM and parallax filtering. Cluster members are clearly separated from the background population, forming a compact over-density.

For the photometric membership filtering of the clusters, we primarily use the G (green) and $G_{\text{RP}}$ (red) passbands from Gaia DR3 (Gaia Collaboration et al. Reference Collaboration2023a). Using these two bands, we created a colour magnitude diagram (CMD) of the stars that passed the astrometric filtering criteria. Upon initial inspection of the CMD, a clear main sequence was identified. We removed a handful of stars on the upper right of the CMD that were several magnitudes from the main sequence and thus very unlikely to be true members even if they were unresolved multiple star systems. There were also a number of stars that appeared slightly below the main sequence or between the white dwarf population and main sequence. We found that most of these outliers were located spatially at the edge of the cluster and towards the edge of the cluster parameter space, with lower membership probabilities. Hence, we manually removed all of these selected stars that did not have parallax and proper motion measurements that were strongly consistent with the rest of the cluster. This ensured there was minimal contamination from non-member stars that passed the initial filtering due to coincidentally having similar astrometric properties to the cluster.

Since we wish to measure the IMF over a broad mass range down to low mass stars, we needed to ensure the completeness of member stars detected by Gaia at faint magnitudes. This was done by comparing identified cluster members from the process described above with those identified using a separate near-infrared survey catalogue. Due to the availability of proper motions, we use observations from the Galactic Clusters Survey (GCS) component of the UKIRT Infrared Deep Sky Survey (UKIDSS) (Lawrence et al. Reference Lawrence2007). Photometric observations were obtained with the Wide Field Camera (Casali et al. Reference Casali2007) using the five infrared filters Z, Y, J, H and K as described in Hodgkin et al. (Reference Hodgkin, Irwin, Hewett and Warren2009). All data was taken from the latest data release UKIDSSDR11PLUS which provides ZYJHK measurements of Praesepe, Pleiades and Alpha Persei with coverages of 28, 79 and 50 sq. degrees respectively.

Figure 2.

Colour magnitude diagrams of filtered Praesepe member stars based on the combined photometric information of Gaia DR3, HIPPARCOS and UKIDSS. Left: Colour magnitude diagram using Gaia G magnitudes and $G-G_{\text{RP}}$ colours. Right: Colour magnitude diagram using UKIDSS J magnitudes and $Y-J$ colours. The colour coding of each star indicates which catalogue it is observed in: either Ga (Gaia), U (UKIDSS), H (HIPPARCOS) or a combination of the three.

When finding cluster member stars from the UKIDSS catalogue, we used a method adapted and slightly changed from Deacon & Hambly (Reference Deacon and Hambly2004) due to the much larger uncertainties in proper motion measurements in the UKIDSS data, compared to Gaia. The method, described in Appendix A of Deacon & Hambly (Reference Deacon and Hambly2004), is a maximum likelihood method which fits cluster stars with a circularly symmetric Gaussian distribution in the proper motion phase space, and field stars with a decaying exponential in the direction of the cluster’s motion and Gaussian perpendicular to it. Initially, stars were filtered using CMDs constructed using photometry in $Z-Y$ and $Y-J$ passbands where the main sequences could be identified using the member stars from the Gaia analysis, cross referenced between catalogues using sky positions. The maximum likelihood was then solved on this subset of stars using a numerical bisection algorithm which can then be converted to give a membership probability of each star as shown in Deacon & Hambly (Reference Deacon and Hambly2004). A star was considered a cluster member if its membership probability was larger than 0.75.

In addition to cross referencing infrared surveys to ensure completeness for low mass stars, we also check the completeness of Gaia for high mass members by comparing observations from its predecessor survey HIPPARCOS (ESA 1997). This serves as a useful comparison as the detectors and optical passbands of HIPPARCOS were able to measure stars with magnitudes as bright as $m_V=-1.44$ without saturating, meaning we can check for and include any bright cluster stars that are beyond the upper limit in brightness of Gaia. The HIPPARCOS catalogue contains proper motions and parallaxes over full sky observations, similar to Gaia, although with much larger uncertainties. We filtered this data similarly to the process described for the Gaia data using the same $\chi^2$ condition. Due to the fewer observations, and thus smaller sample statistics, the mean values used in Equation (1) were found using the HIPPARCOS observations that were matched to the filtered Gaia subset on their sky positions. We then used the photometry in the V and B bands to confirm that all identified cluster members fit along a main sequence in the CMD.

When combining the three catalogues, we prioritise any Gaia measurements and membership probabilities due to their significantly higher accuracy when a star is measured in multiple catalogues. More specifically, if a star was identified as a member in either UKIDSS or HIPPARCOS but the same star was rejected based on its Gaia data, it was still rejected from the final dataset and if a star was accepted based on its Gaia data, it was accepted regardless of the membership results from the other catalogues.

For Praesepe, we find a total of 1 274 member stars. Of these stars, 938 were observed in both Gaia and UKIDSS, 136 only in Gaia, 175 faint stars only in UKIDSS, 23 in both Gaia and HIPPARCOS and 2 stars observed in all three surveys. The magnitudes covered by the different surveys are shown in Figure 2. We observe that Gaia includes the brightest stars in the cluster as well as faint stars all the way down to $Y\approx17$ mag. There are a handful of stars between $Y=12-17$ mag which are only identified in UKIDSS and not Gaia, although there are many more at brighter magnitudes than $G\approx13$ which are only observed in Gaia. This highlights that Gaia is highly complete for the brighter cluster population, whereas UKIDSS provides deeper photometry that extends the detected main sequence to lower masses, particularly at magnitudes where the Gaia sample becomes incomplete.

2.2 Pleiades & Alpha Persei

For the Pleiades we use the same method as for Praesepe to determine cluster members. We again use Gaia DR3 observations within a radius of $15^{\circ}$ around the cluster centre of $\alpha=03^{\text{h}}46^{\text{m}}24^{\text{s}}$ and $\delta=24^{\circ}06'50''$ from Tarricq et al. (Reference Tarricq2021). For the mean proper motions and parallax of the cluster we use starting values of $\bar{\omega}=7.364$ mas, $\bar{\unicode{x03BC}}_\alpha=19.997\,\mathrm{mas\,yr^{-1}}$ and $\bar{\unicode{x03BC}}_\delta=-45.548\,\mathrm{mas\,yr^{-1}}$ from Gaia Collaboration et al. (Reference Collaboration2018). We use the same estimates as Praesepe for dispersion values in the initial Gaia filtering.

After combining catalogues and applying photometric filtering we find a total of 1 602 member stars for the Pleiades. Of these stars, 1 174 are detected in both Gaia and UKIDSS, 226 in UKIDSS only, 148 in Gaia only, 41 in both Gaia and HIPPARCOS, 9 in all three surveys and 4 in HIPPARCOS only.

In the case of Alpha Persei, we slightly modify the filtering process outlined above due to the smaller absolute proper motion of the cluster, which prevents the initial $\chi^2$ filtering from reaching a stable solution for the cluster parameters due to the dense stellar background. The initial dataset was first filtered using the same maximum likelihood method as applied to the UKIDSS data in the previous two clusters. A membership probability threshold of 0.75 was set to remove any very low probability background stars from the Gaia dataset. Since this step removed most of the background stars but still included a significant number of contaminants, we then applied the $\chi^2$ filtering from Equation (1) which was able to reach a stable solution.

Gaia data was again obtained within a $15^{\circ}$ radius of the cluster centre at $\alpha=03^{\text{h}}26^{\text{m}}28^{\text{s}}$ and $\delta=48^{\circ}58'30''$ from Tarricq et al. (Reference Tarricq2021). The initial proper motions and parallax were set to be $\bar{\omega}=5.718$ mas, $\bar{\unicode{x03BC}}_\alpha=19.997\,\mathrm{mas\,yr^{-1}}$ and $\bar{\unicode{x03BC}}_\delta=-45.548\,\mathrm{mas\,yr^{-1}}$ (Gaia Collaboration et al. Reference Collaboration2018).

After merging the catalogues, we found there was a halo of stars around the cluster at radii greater than 6 deg that were only detected in the Gaia dataset. These stars were scattered randomly above and below the apparent main sequence in the cluster CMD and were thus considered non-members and removed. After photometric filtering, we found a total of 1 046 member stars for Alpha Persei. 672 of these members were found in both Gaia and UKIDSS, 201 were only in UKIDSS, 125 were only in Gaia, 45 were in both Gaia and HIPPARCOS and 3 stars were in all three surveys. The raw, uncorrected luminosity functions for all three clusters are provided in Appendix 3 (Table C1 and Figure C1) for reference.

2.3 Summary of astrometric and kinematic parameters

In Table 1, we present the astrometric and kinematic parameters derived in this work from the combined datasets of each cluster. For the cluster centre coordinates $\alpha_0$ and $\delta_0$ , we use the density centre method from von Hoerner (Reference von Hoerner1960, Reference von Hoerner1963) which calculates the average position of stars weighted by the local stellar density. Our derived centres agree well within a few arcminutes to those found by Gaia Collaboration et al. (Reference Collaboration2023b).

For the proper motions $\unicode{x03BC}_{\alpha,{c}}$ and $\unicode{x03BC}_{\delta,{c}}$ , radial velocity $v_{r}$ and parallax $\omega$ , we find the mean and dispersion values of each cluster with uncertainties using the maximum likelihood approach outlined in Pryor & Meylan (Reference Pryor and Meylan1993), which accounts for the individual measurement errors of each star. We also report the mean distance of each cluster d by simply taking the inverse of the mean parallax.

The mean distances we obtain are consistent with existing studies. For the Pleiades, we find a distance of $135.90\pm0.10$ pc confirming newer estimates such as the VLBI result of $136.2\pm1.2$ pc from Melis et al. (Reference Melis, Reid, Mioduszewski, Stauffer and Bower2014), and ruling out the previous HIPPARCOS estimate of $120.2\pm1.9$ pc (Van Leeuwen Reference Van Leeuwen2009). The derived 1D internal velocity dispersions of $\sigma_{v}\approx0.7$ – $0.8\,\mathrm{km\,s^{-1}}$ are typical for nearby open clusters and consistent with earlier findings (e.g. Tarricq et al. Reference Tarricq2021). The mean radial velocities are also in good agreement with Tarricq et al. (Reference Tarricq2021) except for Alpha Persei, which they found to have a low, yet positive value of $0.72\pm0.16\,\mathrm{km\,s^{-1}}$ . From our analysis, we find the cluster to have a slightly negative radial velocity of $-0.76\pm0.27\,\mathrm{km\,s^{-1}}$ , about $1\,\mathrm{km\,s^{-1}}$ different from this previous estimate.

Table 1.

Astrometric and kinematic parameters derived in this work for Praesepe, Pleiades and Alpha Persei. Reported values include the calculated cluster density centre in RA ( $\alpha_0$ ) and DEC ( $\delta_0$ ), mean proper motions in RA ( $\unicode{x03BC}_{\alpha,{c}}$ ) and DEC ( $\unicode{x03BC}_{\delta,{c}}$ ), internal 1D velocity dispersion ( $\sigma_{v}$ ), mean radial velocity ( $v_{r}$ ), mean parallax ( $\omega$ ) with dispersion ( $\sigma_\omega$ ) and distance (d). Uncertainties correspond to 1 $\sigma$ errors.

Table 1. Long description

A table comparing astrometric and kinematic parameters of three star clusters: Praesepe, Pleiades, and Alpha Persei. The table has 8 rows and 10 columns. Column headers are Parameter, Praesepe, Pleiades, and Alpha Persei. Row 1: alpha 0, 08h 40m 12s, 03h 46m 47s, 03h 26m 15s. Row 2: delta 0, +19° 38’ 55’’, +24° 08’ 13’’, +48° 54’ 15’’. Row 3: mu sub alpha cos delta (mas/yr), -31.455 ± 0.034, 12.792 ± 0.019, 18.941 ± 0.030. Row 4: mu sub delta (mas/yr), -11.290 ± 0.028, -29.144 ± 0.023, -21.089 ± 0.031. Row 5: sigma v (km/s), 0.809 ± 0.017, 0.739 ± 0.011, 0.799 ± 0.017. Row 6: v rad (km/s), 34.59 ± 0.19, 5.23 ± 0.26, -0.76 ± 0.27. Row 7: omega (mas), 5.4028 ± 0.0043, 7.3584 ± 0.0055, 5.7188 ± 0.0069. Row 8: sigma omega (mas), 0.110 ± 0.004, 0.177 ± 0.005, 0.175 ± 0.006. Row 9: d (pc), 185.09 ± 0.15, 135.90 ± 0.10, 174.86 ± 0.21.

3.1 Isochrone selection

To estimate the mass of observed cluster member stars, we used PARSEC v2.0 isochrones (Nguyen et al. Reference Nguyen2022) and Mesa Isochrones and Stellar Tracks (MIST) (Choi et al. Reference Choi2016). These models were selected because they incorporate updated input physics, bolometric corrections in the Gaia photometric system, and have parameter grids that extend to the low ages of the observed clusters, including pre-main sequence evolution tracks.

For both models, we adopted previously measured values for the age and metallicity of each cluster. For the Pleiades, we adopt an age of $\approx118\,\mathrm{Myr}$ based on lithium depletion boundary (LDB) measurements, which provide a largely model-independent and precise chronometer for young clusters (e.g. Navascues et al. Reference Navascues, Stauffer and Jayawardhana2004; Frasca et al. Reference Frasca2025). For metallicity, we adopt [M/H] $\approx+0.03$ dex from high resolution spectroscopic analyses of abundances (e.g. Soderblom et al. Reference Soderblom, Laskar, Valenti, Stauffer and Rebull2009; Grilo et al. Reference Grilo2024; Frasca et al. Reference Frasca2025). We also corrected for reddening in the Pleiades isochrone using $A_v = 0.12$ mag (O’Dell et al. Reference O’Dell, Hendry and Collier Cameron1994).

For Alpha Persei, we adopt an age of $\approx80\,\mathrm{Myr}$ , primarily determined by LDB and gyrochronology measurements (e.g. Stauffer et al. Reference Stauffer1999; Navascues et al. Reference Navascues, Stauffer and Jayawardhana2004; Galindo-Guil et al. Reference Galindo-Guil2022; Boyle & Bouma Reference Boyle and Bouma2023). The cluster’s metallicity was taken to be [M/H] $\approx+0.02$ dex based on spectroscopic studies (e.g. Boesgaard et al. Reference Boesgaard, Budge and Ramsay1988; Boyle & Bouma Reference Boyle and Bouma2023). Similar to the Pleiades, we corrected for reddening in Alpha Persei using an extinction of $A_v = 0.30$ mag (Prosser Reference Prosser1992).

For Praesepe, we adopt the parameters derived by Alfonso & García-Varela (Reference Alfonso and García-Varela2023) who report a best-fitting isochrone age of $\log(\mathrm{Age}) = 8.88$ ( $\approx759\,\mathrm{Myr}$ ) and a bulk metallicity of $Z = 0.020$ ([M/H] $\,\approx0.15$ ) from Gaia photometry isochrone fitting. These values are in agreement with other recent isochrone estimates (e.g. Brandt & Huang Reference Brandt and Huang2015; Gaia Collaboration et al. Reference Collaboration2018; Tarricq et al. Reference Tarricq2021).

Figure 3 shows the adopted isochrones overlaid on the cluster CMDs. At intermediate and high stellar masses, the PARSEC and MIST models are in close agreement with one another as well as the observed cluster data. At lower masses, however, we observe that the two models diverge around $M_G\approx7\,\mathrm{mag}$ ( $m\approx0.7\,\mathrm{M}_\odot$ ). These discrepancies likely reflect differences in the models such as the treatment of convection, atmospheric boundary conditions, and low-temperature opacities. The PARSEC models better reproduce the observed photometry of each cluster at low masses, though the variation between the models indicates a degree of uncertainty, particularly for MIST models, in the theory of low mass stellar evolution used in these codes.

Figure 3.

Colour-magnitude diagrams (CMDs) for Alpha Persei, Pleiades, and Praesepe (left to right). The top row shows the absolute magnitude $M_G$ against ( $G-G_{\text{RP}}$ ), while the bottom row shows $M_G$ against ( $G_{\text{BP}}-G_{\text{RP}}$ ). In each panel, the PARSEC v2.0 isochrone using values from prior work is shown in orange and the corresponding MIST isochrone in blue. Grey dashed lines indicate masses along the isochrone for PARSEC models. Some The comparison illustrates the agreement between the two stellar evolution models at higher masses but some differences in photometry for stars fainter than $M_G\approx7$ mag.

Stellar masses were assigned via isochrone interpolation in three photometric bands: $G, G_{\text{BP}}$ , and $G_{\text{RP}}$ . For each star, we computed a $\chi^2$ statistic of the form

(3) \begin{equation}\chi^2(m) = \sum_{i \in \{G, G_{\rm BP}, G_{\rm RP}\}}\frac{(M_{i,\mathrm{obs}} - M_{i,\mathrm{iso}}(m))^2}{\sigma_{i}^2},\end{equation}

where $M_{i,\mathrm{obs}}$ and $\sigma_i$ are the observed magnitudes and uncertainties, and $M_{i,\mathrm{iso}}(m)$ is the model magnitude at mass m along the isochrone. The $\chi^2$ values were converted into relative likelihoods via

(4) \begin{equation}\mathcal{L}(m) \propto \exp\left(-\frac{\chi^2(m)}{2}\right),\end{equation}

which were then normalised to obtain a posterior probability distribution function (PDF) in mass under the assumption of uniform priors. The adopted stellar mass corresponds to the median value of this PDF, and the $1\sigma$ uncertainty was defined from the central 68% interval.

This procedure was repeated independently for the PARSEC and MIST grids, enabling a direct comparison of derived masses and an assessment of model-dependent uncertainties in the resulting PDMF.

3.2 Empirical mass-luminosity relation using eclipsing binaries

As there is evidence that the theoretical mass conversion from the two investigated isochrone models, particularly for low mass stars, may be inconsistent or uncertain, we choose to estimate stellar masses using an independent, empirically determined mass-luminosity (M-L) relationship. To achieve this, we use a similar method to Kroupa et al. (Reference Kroupa, Tout and Gilmore1993) who fit a model to mass and magnitude measurements of eclipsing binary systems from Popper (Reference Popper1980). In this work, we use a collection of recent eclipsing binary measurements with a focus on low mass stars since we are more confident in the isochrone modelling for high mass stars $\gt 1\,\mathrm{M}_\odot$ . The individual star data used is tabulated in Table A1; the majority of these observations come from previous collections compiled in either Torres et al. (Reference Torres, Andersen and Giménez2010), Benedict et al. (Reference Benedict2016) or Iglesias-Marzoa et al. (Reference Iglesias-Marzoa, López-Morales, Arévalo, Coughlin and Lázaro2017).

It is important to note that this relationship is determined using $M_V$ measurements as it has historically been the standard optical passband for eclipsing binary analyses and nearby benchmark stars, and is not directly equivalent to Gaia G-band photometry. Thus, we first compute the fit between mass and $M_V$ using the binary star data and then convert the Gaia observations of each cluster to $M_V$ using the colour-dependent transformation provided by Riello et al. (Reference Riello2021), based on the Gaia EDR3 photometric calibration.

The comprehensive catalogue presented in Torres et al. (Reference Torres, Andersen and Giménez2010) extends over a large range of stellar masses. Since we are only interested in main sequence stars, particularly with low masses, we exclude stars brighter than $M_V=3$ ( $m\approx1.4 \,\mathrm{M}_\odot$ ) which have undergone post main-sequence evolution. We also exclude the 9 stars from this dataset with magnitudes fainter than $M_V=7.5$ ( $m\approx0.7\,\mathrm{M}_\odot$ ) as the adopted polynomial function of bolometric corrections (BC) from Flower (Reference Flower1996) begins to break down at very low masses.

To better populate the parameter space with stars between $0.3$ – $0.5\,\mathrm{M}_\odot$ we apply our own bolometric corrections to the bolometric magnitudes of six stars from Birkby et al. (Reference Birkby2012). To achieve this, we used observations of O9-M9 dwarf effective temperatures and BC values as determined in Pecaut & Mamajek (Reference Pecaut and Mamajek2013) and used them to interpolate BC values based on the observed $T_{\rm eff}$ values in Birkby et al. (Reference Birkby2012), with uncertainties being propagated to the final $M_V$ value.

With the data compiled, we fitted a smoothing B-spline satisfying the Generalised Cross Validation (GCV) criterion as described in Wahba (Reference Wahba1990). This relationship with respect to the eclipsing binary data and the previous fit from Kroupa et al. (Reference Kroupa, Tout and Gilmore1993) are displayed in Figure 4. We see similar behaviour between the two models, although our model does appear to fit our sample better between around $0.5$ – $0.9\,\mathrm{M}_\odot$ and flattens more at magnitudes fainter than $M_V\approx16$ .

Figure 4.

Empirical mass–luminosity relation for low-mass stars. The black points show the sample of eclipsing binaries used to calibrate the relation given in Table A1. The blue curve shows the spline fit derived in this work, while the dashed red line shows the spline fit from Kroupa et al. (Reference Kroupa, Tout and Gilmore1993). The spline fit provides a smooth empirical calibration of the mass–luminosity relation across the sampled magnitude range and is adopted in this work to estimate stellar masses from $M_V$ .

While using an empirical M-L relation avoids issues related to theoretical modelling of stellar evolution, it does have several limitations in the context of star cluster mass estimates. Unlike the multi-band fitting described above, which simultaneously uses $G, G_{\text{BP}}$ , and $G_{\text{RP}}$ measurements with errors, the $M_V$ relation does not account for colour information and is therefore less sensitive to subtle temperature or metallicity effects. In fact, of the low mass stars with metallicity estimates collated in Benedict et al. (Reference Benedict2016), we find a moderately positive relationship between the residuals of our model fit and stellar metallicity with a Pearson correlation coefficient of $r=0.51$ that is statistically significant ( $p=0.0023$ ). This would imply a star with higher metallicity may have a larger mass than one with lower metallicity and the same $M_V$ . We do not, however, incorporate corrections for metallicity in this work as robust metallicity estimates are unavailable for the majority of stars in our binary sample. The potential systematic bias introduced by neglecting metallicity dependence should therefore be borne in mind when interpreting the derived mass function.

3.3 Mass function modelling and fitting

We infer the stellar mass function (MF) using a forward modelling approach that explicitly accounts for measurement uncertainties on individual stellar masses. We choose to model the mass function as a multistage power law function similar to Kroupa (Reference Kroupa2001), for better direct comparisons with determinations of IMF values from other published work. This function is commonly broken into either two or three stages depending on the particular mass range being fitted. A three stage power law MF can be described by the function

(5) \begin{equation} \xi(m) = \begin{cases} Cm^{-\alpha_{l}}, & m_{\text{min}} \leq m \leq m_{x1}\\ Cm_{x1}^{\alpha_{m} - \alpha_{l}}m^{-\alpha_{m}}, & m_{x1} \leq m \leq m_{x2}\\ C m_{x2}^{\alpha_{h} - \alpha_{m}} m_{x1}^{\alpha_{m} - \alpha_{l}} m^{-\alpha_{h}}, & m_{x2} \leq m \leq m_{\text{max}} \end{cases} \end{equation}

where C is a normalisation constant, $m_{x1}$ and $m_{x2}$ are the break masses between the different MF segments, $m_{\text{min}}$ and $m_{\text{max}}$ are the minimum and maximum masses, and $\alpha_{i}$ is the mass function ‘slope’ of mass domain i. Here, the constants $m_{x1}^{\alpha_{m} - \alpha_{l}}$ and $m_{x2}^{\alpha_{h} - \alpha_{m}} m_{x1}^{\alpha_{m} - \alpha_{l}}$ ensure continuity between the different power law segments. A two stage power law MF is similar in form, with only the first two segments, needing only two slope values and a single break point to be defined.

In general, the parameter of interest is the MF slope $\alpha_i$ as it describes the relative amount of high and low mass stars in a given mass range i. Here we describe the fitting process in which we determined the best fitting two and three stage power law functions for each cluster in a Bayesian framework, and find the preferred model given the Bayes factor K. We applied this method to the stellar masses estimated using both PARSEC and MIST isochrone models as well as the empirical M-L spline fit determined in Section 3.2 and compare the various forms and values of the PDMF inferred using each model.

If $m_j^{\text{obs}}$ and $\sigma_j$ denote the observed mass and uncertainty of star j, then instead of treating each mass as an exact sample from a MF $\xi(m)$ , we marginalise over the unobserved true mass $m_j^{\text{true}}$ . The likelihood contribution of each observed star is then

(6) \begin{equation} \mathcal{L}_j = \int_{m_{\text{min}}}^{m_{\text{max}}}\xi(m)P \left(m_j^{\text{obs}} | m,\sigma_j \right)dm,\end{equation}

where $P(m_j^{\text{obs}} | m,\sigma_j)$ is assumed to be a Gaussian probability density centred at m with dispersion $\sigma_j$ .

Figure 5.

Posterior distributions for the parameters of the three-segment broken power-law IMF model for the Pleiades. The fitted parameters include the three power-law slopes ( $\alpha_1$ , $\alpha_2$ , $\alpha_3$ ), the two break masses ( $m_{\mathrm{break},1}$ and $m_{\mathrm{break},2}$ ), and the minimum and maximum stellar masses ( $m_{\min}$ and $m_{\max}$ ) with $1\sigma$ credible intervals. Most parameters, excluding $m_{\max}$ , are well constrained with approximately uni-modal posterior distributions, while several moderate correlations are evident between the power-law slopes and break masses.

We compute this integral numerically over the intersection of the model bounds $\left[m_{\text{min}}, m_{\text{max}}\right]$ and a finite window $\big[m_j^{\text{obs}}\pm5\sigma_j\big]$ using a trapezoidal method on a fixed grid of points. The total log-likelihood is then simply given by

(7) \begin{equation} \ln\mathcal{L}=\sum_{j=1}^N\ln\mathcal{L}_j.\end{equation}

This leaves us with 5 free parameters for a two stage power law fit ( $\alpha_1$ , $\alpha_2$ , $m_x$ , $m_{\text{min}}$ , $m_{\text{max}}$ ) and 7 free parameters for a three stage power law fit ( $\alpha_1$ , $\alpha_2$ , $\alpha_3$ , $m_{x1}$ , $m_{x2}$ , $m_{\text{min}}$ , $m_{\text{max}}$ ) that need to be determined. We obtain posterior distributions and the Bayesian evidence of each model using the nested sampling algorithm implemented in NAUTILUS (Lange Reference Lange2023), assuming uniform priors for each parameter. An example of a posterior distribution obtained using this method is given in Figure 5, where we simultaneously fit all parameters of a three stage power law to the Pleiades masses found using the PARSEC isochrone model. We see good convergence of all parameters with approximately Gaussian errors, excluding the maximum mass bound $m_{\text{max}}$ .

Figure 6.

Present-day mass functions (PDMFs) for Alpha Persei, the Pleiades, and Praesepe (left to right). In each panel, the favoured broken power-law model (either two- or three-segment, depending on the cluster) is shown for masses estimated using three different methods: PARSEC isochrones (solid blue), MIST isochrones (dashed orange), and the empirical mass–luminosity relation (dot–dashed green). The resulting PDMF shapes are broadly consistent across the three mass estimation methods, however, the models display differences below $\approx0.3\,\mathrm{M}_\odot$ . For Praesepe (right) there was no significant preference between models using PARSEC masses so both the three segment (blue line) and two segment (blue dotted) are shown.

We apply a manual lower bound to each cluster corresponding to ${G}\approx20\,\mathrm{mag}$ for Alpha Persei and Pleiades and ${G}\approx19.2\,\mathrm{mag}$ for Praesepe due to the completeness of Gaia at fainter magnitudes as well as the limitations in the isochrones at any smaller masses. This corresponds to a lower limit of $0.10\,\mathrm{M}_\odot$ for Alpha Persei and Pleiades and $0.17\,\mathrm{M}_\odot$ for Praesepe.

When determining the optimal functional form for the PDMF, we find varying results depending on the mass estimation model used. For example, in the Pleiades we find very strong evidence ( $\log K\gg2$ ) for a three stage fit when using either PARSEC or MIST mass estimates but find strong evidence ( $\log K\gt 1$ ) for a two stage fit when using the empirical M-L relation mass estimates. We find a similar trend between the three mass estimates for Alpha Persei, with the two isochrone models producing similar three stage fits and the spline fit model giving preference to a two stage fit, similar to that of the Pleiades. Praesepe, however, appears to be unique to the other clusters in that the MIST isochrone masses and the empirical M-L give evidence for a two stage PDMF with similar slope and break values. The PARSEC models, on the other hand, give no preference ( $\log K\lt 0.2$ ) to a 2 stage or 3 stage fit with slightly different parameter values compared to the other two fits. In this case, we choose to use the 3 stage PARSEC fit for Praesepe so that it more easily comparable to values of the other two clusters.

The differences between the three mass estimates for each cluster can be seen in Figure 6. We see that the MF slope and break mass at high masses $\approx1\,\mathrm{M}_\odot$ is consistent between the three models, giving confidence in our estimates in that mass range. This agreement begins to break down at masses $\lesssim0.5\,\mathrm{M}_\odot$ , where the MF from different mass estimation models begin to diverge. Nevertheless, the broad structure of the PDMF, with a steep high-mass slope, a flattening toward lower masses, and a break near $m_x\approx1\,\mathrm{M}_\odot$ is recovered by all three approaches. The empirical two-stage fits (Table 2) show remarkably consistent break masses across all clusters with $m_x\approx0.86$ – $0.88\,\mathrm{M}_\odot$ and low-mass slopes in the range $\alpha_\mathrm{low}\approx1.2$ – $1.4$ , while the high-mass slopes vary at a $\approx0.3$ – $0.5$ level. These similarities suggest that the global PDMF shape is generally robust, whereas the detailed behaviour at lower masses is sensitive to the adopted mass–luminosity prescription.

Given these considerations, we chose to use the PARSEC mass estimates for the rest of the analysis of each cluster. The PARSEC models provide three stage solutions for all three clusters with statistical evidence, allowing us to be consistent in the parametrisation of the MF. Furthermore, we note that PARSEC produces better CMD fits compared to the MIST isochrones, incorporates pre-main sequence evolution and considers the age and metallicity dependence of a star’s colour and luminosity. The empirical relationship, however, averages over stars with a range of metallicities and ages and does not consider pre-main sequence evolution which could be important for the low mass stars in Alpha Persei and Pleiades in particular. We therefore use PARSEC as our mass estimation, while acknowledging that systematic differences at $m \lesssim 0.5\,\mathrm{M}_\odot$ introduce model-dependent variations in the detailed low-mass behaviour. The best fitting PDMF models for PARSEC derived masses can be seen compared to the data in Figure 7 with parameter estimations given in Table 2.

Table 2.

Best-fitting present-day mass function (PDMF) parameters of Alpha Persei, Pleiades, and Praesepe using PARSEC v2.0 isochrone and empirical mass-luminosity relation mass estimates. Reported are the observed three-stage uncorrected PDMF slopes using PARSEC ( $\alpha_{h,{\rm pc}}$ , $\alpha_{m,{\rm pc}}$ , $\alpha_{l,{\rm pc}}$ ) with mass break points ( $m_{x1, {\rm pc}}$ , $m_{x2, {\rm pc}}$ ) as well as the two-stage PDMF slopes using the empirical M-L masses ( $\alpha_{h,{\rm ML}}$ , $\alpha_{l,{\rm ML}}$ ) with break mass $m_{x,{\rm ML}}$ .

Table 2. Long description

A table comparing mass function parameters for Alpha Persei, Pleiades, and Praesepe. The table has 8 rows and 4 columns. The columns are labeled Parameter, Alpha Persei, Pleiades, and Praesepe. The rows list different parameters and their corresponding values for each cluster. Row 1: Parameter, Alpha Persei, Pleiades, Praesepe. Row 2: αh,pc, 2.28 ± 0.38, 3.04 ± 0.34, 2.53 ± 0.85. Row 3: αm,pc, 1.73 ± 0.24, 1.62 ± 0.13, 1.99 ± 0.47. Row 4: αl,pc, 0.32 ± 0.34, -0.29 ± 0.30, 0.84 ± 0.31. Row 5: mh,pc (M☉), 0.95 ± 0.34, 0.91 ± 0.12, 1.20 ± 0.33. Row 6: ml,pc (M☉), 0.30 ± 0.06, 0.24 ± 0.02, 0.50 ± 0.11. Row 7: αh,ML, 2.52 ± 0.20, 3.19 ± 0.24, 2.80 ± 0.40. Row 8: αML, 1.29 ± 0.08, 1.36 ± 0.06, 1.22 ± 0.08. Row 9: mh,ML (M☉), 0.86 ± 0.̄08, 0.88 ± 0.07, 0.88 ± 0.11.

Figure 7.

Cluster observed present day mass functions using PARSEC v2.0 isochrone masses, uncorrected for binaries or dynamical evolution. Black points show the fraction of stars in each log mass bin with $1\sigma$ error bars, assuming Poisson statistics. The solid red line indicates the best-fitting three-segment broken power-law model.

We quote the observed total mass $M_\mathrm{tot}$ by simply summing all the estimated masses for each member star, and average mass $\bar{M}$ by taking the mean. The projected 2D half mass radius $r_{h,{\rm 2D}}$ was found by first finding each star’s radius from the cluster density centre given in Table 1 and then sorting the stars in order of radius to construct the cumulative mass distribution. We identify $r_{h,{\rm 2D}}$ as the radius at which the cumulative is half the total mass, using linear interpolation between the two nearest points if the half mass is between two cumulative mass points. For the 3D half mass radius, we first found the heliocentric Cartesian coordinates of each member star using their right ascension, declination and parallax. We then determined the three dimensional cluster centre, weighted by stellar mass and found the radius of each cluster member. The half mass radius was then determined similarly to the 2D projected case. We estimated uncertainty values for $M_\mathrm{tot}$ , $\bar{M}$ , $r_{h,{\rm 2D}}$ and $r_{h,{\rm 3D}}$ by resampling from all cluster members and bootstrapping over 50 000 iterations. The standard deviation of these measurements over all iterations describe the sensitivity of these measurements to member selection and thus are used as the quoted errors.

The half-mass relaxation time $t_\mathrm{rh}$ of each cluster was estimated using the derivation from Spitzer (Reference Spitzer1987), which is given by

(8) \begin{equation} t_\mathrm{rh} = 0.138\frac{M_\mathrm{tot}^{1/2}r_{h,{\rm 3D}}^{3/2}}{G^{1/2}\bar{M}\ln\left(\gamma\frac{M_\mathrm{tot}}{\bar{M}}\right)} \text{Myr}\end{equation}

where G is the gravitational constant and $\gamma=0.11$ (Giersz & Heggie Reference Giersz and Heggie1996). Uncertainties from $M_\mathrm{tot}$ , $\bar{M}$ and $r_{h,{\rm 3D}}$ were also propagated to give an error estimate for each cluster. Due to its small age and large observed half mass radius, Alpha Persei is younger than one relaxation time, indicating that it has likely not been strongly affected by dynamical evolution. Pleiades and Praesepe on the other hand have ages similar to (Pleiades) or much greater (Praesepe) than their relaxation times, meaning that correcting for dynamical effects is more important when estimating the initial states of these clusters. A summary of the derived values based on estimated stellar masses is presented in Table 3.

Table 3.

Observed structural parameters of Alpha Persei, Pleiades, and Praesepe. Reported values include the total cluster masses ( $M_\mathrm{tot}$ ), mean stellar masses ( $\bar{M}$ ), projected 2D half-mass radii ( $r_{h,{\rm 2D}}$ ), 3D half-mass radius ( $r_{h,{\rm 3D}}$ ) and relaxation times ( $t_\mathrm{rh}$ ). Uncertainties correspond to 1 $\sigma$ errors.

Table 3. Long description

A table comparing structural parameters of Alpha Persei, Pleiades, and Praesepe. The table has 4 rows and 5 columns. The columns are labeled Value, Alpha Persei, Pleiades, and Praesepe. The rows are labeled with different parameters: M_tot (M_sun), M̄ (M_sun), r_h,2D (pc), r_h,3D (pc), and t_rh (Myr). Each cell contains numerical values with uncertainties. Row 1: M_tot (M_sun), Alpha Persei 518 ± 21, Pleiades 650 ± 18, Praesepe 598 ± 14. Row 2: M̄ (M_sun), Alpha Persei 0.62 ± 0.03, Pleiades 0.50 ± 0.02, Praesepe 0.574 ± 0.02. Row 3: r_h,2D (pc), Alpha Persei 5.93 ± 0.33, Pleiades 3.02 ± 0.14, Praesepe 3.16 ± 0.15. Row 4: r_h,3D (pc), Alpha Persei 6.89 ± 0.29, Pleiades 4.27 ± 0.13, Praesepe 4.92 ± 0.17. Row 5: t_rh (Myr), Alpha Persei 305 ± 24, Pleiades 190 ± 11, Praesepe 208 ± 12.

4.1 Observable binaries

For the observed clusters, we measured the stick-out binary fraction by first modelling the main sequence below the turnoff as a high order polynomial to better match the observations of low mass stars compared to the theoretical isochrone. We then identified systems that deviated from the main sequence by 2.5 $\sigma_{p}$ or more to determine the observed binary fraction of each cluster, where $\sigma_{p}$ is a star’s photometric error in colour-magnitude space. When quoting any binary fractions, we refer to the fraction defined by

(9) \begin{equation} f_\mathrm{bin} = \frac{N_\mathrm{bin}}{N_\mathrm{sing} + N_\mathrm{bin}}\end{equation}

where $N_\mathrm{bin}$ is the number of binary systems and $N_\mathrm{sing}$ is the number of single stars.

Figure 8 depicts the identified stick-out binaries in Praesepe. These systems are more concentrated toward the cluster centre, as expected for genuine binary members. Such a distribution is consistent with the effects of dynamical mass segregation, where these more massive systems would migrate inward over time. This spatial distribution therefore supports the interpretation that these stick-out stars are indeed binary systems which are cluster members, rather than contaminating stars with proper motions that coincidentally align with that of the cluster.

Figure 8.

Identified ‘stick-out’ stars in Praesepe. (left) CMD of Praesepe with observable ‘stick-out’ stars that lie above the main sequence highlighted with red triangles and the main sequence polynomial fit in blue. The CMD uses Gaia photometry in the G and $G_{\text{RP}}$ passbands. (centre) Praesepe stellar spatial distribution with stick-out stars highlighted in red. (right) Cumulative fraction of stick-out stars (red dashed) and main sequence stars (black) as a function of radius. The stick-out stars are strongly concentrated to the cluster centre, indicating they are likely cluster members rather than contaminants.

Using this method, we observe the stick-out binary fraction, defined as $f_{\mathrm{SO}}=N_{\mathrm{SO}}/(N_\mathrm{MS}+N_{\mathrm{SO}})$ , of Alpha Persei, Pleiades and Praesepe to be $(18.3\pm1.4)\%$ , $(16.4\pm1.2)\%$ and $(18.8\pm1.3)\%$ respectively. The uncertainties quoted here were estimated using bootstrapping and resampling from cluster members with replacement.

4.2 Synthetic cluster creation

To recover the true binary fraction of the cluster, we must also consider binary systems that will still lie close to the main sequence due to one star being much more massive and brighter than the other star in a given pair. We did this by creating synthetic stellar populations and CMDs in a similar procedure as in Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013), however we implement the method in a Bayesian framework rather than a direct grid search. Thus we predict the binary corrected MF and present-day binary fraction of each cluster using the following procedure.

We first drew $N_\mathrm{sys}$ masses from a given MF, defined by three mass function slopes $\alpha_{l}$ , $\alpha_{m}$ and $\alpha_{h}$ and best fitting break points $m_{x1}$ and $m_{x2}$ as in Table 2 for each cluster using inverse sampling, where random uniform values are mapped through the cumulative distribution function of the MF to generate masses. The photometry of each star was then derived using the PARSEC isochrones and interpolating $G, G_\mathrm{BP}$ and $G_\mathrm{RP}$ magnitudes for each mass. Of these stars, we selected $N_\mathrm{bin}$ stars to be the primary components of binaries, where $N_\mathrm{bin}=N_{\text{star obs}}\cdot f_\mathrm{bin}$ . For each primary component, a secondary mass was selected from a flat mass ratio distribution between $m_\mathrm{min}$ and $m_\mathrm{primary}$ . The photometry of each secondary mass was derived similarly using the isochrone models and combined with the magnitudes of the primary components to give the total magnitude of each binary system as

(10) \begin{equation} {G}_\mathrm{bin} = -2.5\log(10^{-0.4{G}_\mathrm{prim}}+10^{-0.4{G}_\mathrm{sec}}),\end{equation}

(11) \begin{equation} {G}_{\text{BP, bin}} = -2.5\log(10^{-0.4{G}_{\text{BP, prim}}}+10^{-0.4{G}_{\text{BP, sec}}}),\end{equation}

(12) \begin{equation} {G}_{\text{RP, bin}} = -2.5\log(10^{-0.4{G}_{\text{RP, prim}}}+10^{-0.4{G}_{\text{RP, sec}}}).\end{equation}

Synthetic photometric errors were introduced using a Gaussian distribution with $\sigma_\mathrm{mag}$ equal to the median error of the photometry from the observed sample of stars.

With this synthetic CMD, we then identified the fraction of stars that deviate significantly from the theoretical isochrone of the cluster, labelling them as stick-out binaries. We also measured the mass function slopes of the synthetic cluster with unresolved binaries so that it could be compared to the observed clusters. For a given set of tested input parameters, we created and measured 5 synthetic clusters to account for the variability in cluster generation.

This synthetic cluster creation method was then implemented in an algorithm to estimate the Bayesian posterior and evidence and find the total binary fraction and mass function that best recreate observations of each cluster. We make use of the nested sampling algorithm in the python package NAUTILUS (Lange Reference Lange2023). Using this approach, we aim to find the combination of PDMF slopes $\alpha_{l, {\rm bin}}$ , $\alpha_{m, {\rm bin}}$ , $\alpha_{h, {\rm bin}}$ and binary fraction $f_\mathrm{bin}$ that maximise the log likelihood given by

(13) \begin{align} \begin{split} \ln\mathcal{L} = &-0.5 \sum_i \Bigg[ \frac{(f_{\text{SO}, i} - f_{\mathrm{SO,obs}})^2}{\sigma_{\mathrm{SO,obs}}^2} + \frac{(\alpha_{{l}, i} - \alpha_{l,\mathrm{obs}})^2}{\sigma_{\alpha_{l}, i}^2 + \sigma_{\alpha_{l},\mathrm{obs}}^2} \\ &+ \frac{(\alpha_{{m}, i} - \alpha_{m,\mathrm{obs}})^2}{\sigma_{\alpha_{m}, i}^2 + \sigma_{\alpha_{m},\mathrm{obs}}^2} + \frac{(\alpha_{{h}, i} - \alpha_{h,\mathrm{obs}})^2}{\sigma_{\alpha_{h}, i}^2 + \sigma_{\alpha_{h},\mathrm{obs}}^2} \Bigg] \end{split}\end{align}

where i is the simulation trial number (1–5 here), $f_{\text{SO}, i}$ is the stick-out fraction of trial i, $f_{\mathrm{SO,obs}}$ is the observed cluster stick-out fraction with uncertainty $\sigma_{\mathrm{SO,obs}}$ , $\alpha_{{x}, i}$ is the apparent mass function slope of trial i with uncertainty $\sigma_{\alpha_{x}, i}$ in the mass range $x\epsilon\{{l,m,h}\}$ when synthetic binaries are unresolved and $\alpha_{x,\mathrm{obs}}$ is the observed cluster mass function with uncertainty $\sigma_{\alpha_{x},\mathrm{obs}}$ in mass range x.

For computational efficiency, we measure the simulated MF in each iteration of the sampler using the maximum likelihood estimator as described in Appendix A1 from Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013), who show that the maximum likelihood estimation (MLE) $\hat{\alpha}$ for a set of masses $x_i$ bounded by $(x_\mathrm{min}, x_\mathrm{max})$ is given by the equation

(14) \begin{equation} \hat{\alpha} = 1+n\left(\sum_{i=1}^n\ln\frac{x_i}{x_\mathrm{min}}-n\frac{\ln X}{1-X^{\hat{\alpha}-1}}\right)^{-1} \end{equation}

where n is the number of stars and $X=\frac{x_\mathrm{max}}{x_\mathrm{min}}$ .

Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013) also show that the error of $\hat{\alpha}$ can be estimated using the negative of the inverse of the Fisher information which gives

(15) \begin{equation} \sigma(\hat{\alpha})=\frac{1}{\sqrt{n}}\left[(\hat{\alpha}-1)^{-2}-\ln^2X\frac{X^{\hat{\alpha}-1}}{(1-X^{\hat{\alpha}-1})^2}\right]^{-1/2}.\end{equation}

Here it is important to make the distinction between the stellar MF and system MF, which are related yet represent different quantities. The system MF refers to the mass distribution of single stars and primary stars of multiple systems, whereas the stellar MF describes the mass distribution of all stars in the population. In this method, we use the system MF to generate the population of single and primary mass stars for binary systems. We then introduce additional low and medium mass stars as secondary stars from flat distributions rather than the same system MF. This results in the binary corrected stellar MF, which describes the distribution of all individual star masses in the cluster. While the system MF is useful in the creation of stellar populations when the properties of binaries are of interest, the stellar MF is generally the preferred description as it fully measures the complete mass distribution. Thus, we convert the best fitting system PDMF into a stellar PDMF by generating a synthetic cluster with a large number of stars to reduce statistical uncertainty, using the system PDMF and binary fraction and then measure the stellar PDMF using Equation (14).

The adopted binary pairing method provides a practical framework that preserves simple analytic forms for both the stellar IMF and the binary mass-ratio distribution. It should be noted that drawing secondary masses from a flat mass-ratio distribution introduces a small excess of systems near $m_\mathrm{min}$ at high binary fractions. However, we find that even for star counts of $N_\mathrm{star}\approx10\,000$ with 100% binary fractions, that the small deviation from a perfect power law is dominated by the scatter from stochastic sampling and the resulting stellar mass function retains its overall three-stage power law form with the same break masses. Quantitatively, the introduced deviations are smaller than the stochastic scatter expected from finite sampling, with KS tests indicating that synthetic populations remain consistent with a three-stage power law description. While alternative pairing schemes are possible (e.g. Kouwenhoven et al. Reference Kouwenhoven, Brown, Goodwin, Portegies Zwart and Kaper2009; Belloni et al. Reference Belloni, Askar, Giersz, Kroupa and Rocha-Pinto2017), the adopted approach provides a simple, computationally efficient framework that can be described analytically and adequately reproduces the stellar populations modelled in this work.

Figure 9 shows the posterior distribution of each cluster. We see that there is good convergence for the estimates of the system MF and binary fraction. The best fitting parameters given by these distributions are listed in Table 4. We find that the calculated binary fractions are compatible with being the same for all three clusters. Furthermore, the intermediate mass and high mass MF slopes are also compatible, with the low mass slopes significantly deviating from one another.

Figure 9.

Posterior distributions for the present day binary fraction and three binary corrected system PDMF slopes for Alpha Persei, Pleiades, and Praesepe. All three clusters have comparable binary fractions $f_\mathrm{bin}$ , high MF slopes $\alpha_{h,{\rm bin}}$ and intermediate MF slopes $\alpha_{m,{\rm bin}}$ but have distinct low MF slopes $\alpha_{l,{\rm bin}}$ .

The binary fractions inferred here for Alpha Persei, the Pleiades, and Praesepe are broadly consistent with values reported in existing studies. For Alpha Persei, we measure $f_\mathrm{bin}=0.238\pm0.012$ . Using similar Monte Carlo modelling, Sheikhi et al. (Reference Sheikhi2016) found a binary fraction of $0.23\pm0.09$ across the mass range $0.1$ – $4.7\,\mathrm{M}_\odot$ , which is in good agreement with our result. Our estimate for the Pleiades, $f_\mathrm{bin}=0.200\pm0.008$ , lies between reported values in previous works. For example, the photometric analysis of Liu et al. (Reference Liu, Shao and Li2025) finds a higher binary fraction of $0.34\pm0.02$ , while Cordoni et al. (Reference Cordoni2023) report a lower value of $0.148\pm0.100$ . Given the different methodologies and sensitivities to binary mass ratios in these studies, our result falls comfortably within the range of existing estimates. Finally, for Praesepe we obtain $f_\mathrm{bin}=0.221\pm0.010$ , which is below some prior estimates such as the $0.35\pm0.05$ binary fraction inferred in Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013), but more consistent with other photometric studies such as the fraction of $0.251\pm0.108$ determined in Cordoni et al. (Reference Cordoni2023). The differences in estimations likely arise due to the various observational methods, stellar mass ranges and detectable mass-ratio limits across different studies.

The effects of binaries on the Pleiades mass function are visually seen in Figure 10. Factoring in unresolved binaries is shown to flatten the mass function slope in the low mass range and slightly steepen the slope for intermediate mass stars with little effect at the high mass end. This behaviour is expected, and the fact that we see the same effect helps confirm the validity of the results using this Bayesian approach. We see only slight variations in the intermediate mass range, meaning the measurement of the MF for stellar samples between $\gtrsim0.25\,\mathrm{M}_\odot$ are not significantly affected by binary or multiple star contaminants. The largest effect is on the MF below $\approx0.25$ , which already has additional uncertainty due to the dependence on the mass estimation model used.

Table 4.

Binary corrected cluster IMF properties. Shown are the estimated binary fraction ( $f_\mathrm{bin}$ ) and the high, intermediate, and low mass PDMF slopes ( $\alpha_{h, {\rm bin}}$ , $\alpha_{m,{\rm bin}}$ , $\alpha_{l,{\rm bin}}$ ) after correcting for unresolved binaries. Reported uncertainties correspond to 1 $\sigma$ errors.

Table 4. Long description

A table comparing binary corrected cluster IMF properties for Alpha Persei, Pleiades, and Praesepe. The table has 4 rows and 4 columns. The columns are labeled ’Value’, ’Alpha Persei’, ’Pleiades’, and ’Praesepe’. The rows are labeled with different values: f_bin, alpha_m, bin, alpha_m, int, and alpha_m, low. Each cell contains numerical values with uncertainties. Row 1: f_bin, 23.8 plus or minus 1.2 percent, 20.0 plus or minus 0.8 percent, 22.1 plus or minus 1.0 percent. Row 2: alpha_m, bin, 2.37 plus or minus 0.20, 3.11 plus or minus 0.19, 2.66 plus or minus 0.45. Row 3: alpha_m, int, 1.76 plus or minus 0.21, 1.69 plus or minus 0.08, 2.01 plus or minus 0.24. Row 4: alpha_m, low, 0.90 plus or minus 0.19, 0.46 plus or minus 0.18, 1.41 plus or minus 0.16.

Figure 10.

Binary corrections of the PDMF of Pleiades. Binary correction leads to a relatively larger fraction of low and intermediate mass stars compared to high mass stars due to additional, unseen low mass stars that are binary companions.

5.1 Emulator training

We ran 942, 672 and 550 complete N-body simulations for Alpha Persei, Pleiades, and Praesepe, respectively, using NBODY7 (Nitadori & Aarseth Reference Nitadori and Aarseth2012). The initial positions and velocities of the simulated clusters were found using the observed central positions, proper motions and radial velocities listed in Table 1, and tracing back the path of the cluster in the galactic potential from Irrgang et al. (Reference Irrgang, Wilcox, Tucker and Schiefelbein2013) using a Leapfrog integrator in GALPY (Bovy Reference Bovy2015). The clusters were generated with a King (Reference King1962) density profile which can be described using a 3D half mass radius $r_{h,{\rm ini}}$ and a dimensionless concentration parameter c. For a King model, the tidal radius $r_t$ is related to the core radius $r_{c}$ through

(16) \begin{equation} r_t=r_{c}10^c.\end{equation}

The secondary masses for primordial binaries were sampled from a flat distribution of mass ratios, similar to the sampling described in Section 4. Semi-major axes were sampled from a log-uniform distribution bounded between a minimum separation of 50 times the combined stellar radius and up to 1 000 AU. These limits would therefore exclude any very close binary systems and wide binary pairs from the initial cluster state in favour of numerical stability and efficiency. Thus, the true initial fraction of binaries may be larger than that explicitly represented within the adopted separation range, since the imposed limits preferentially exclude the closest and widest systems that contribute to the full observationally inferred binary population. Alternative prescriptions for constructing primordial binary populations, including different pairing algorithms and orbital parameter distributions motivated by observations and star formation models, have also been proposed (e.g. Kouwenhoven et al. Reference Kouwenhoven, Brown, Goodwin, Portegies Zwart and Kaper2009; Moe & Di Stefano Reference Moe and Di Stefano2017; Belloni et al. Reference Belloni, Askar, Giersz, Kroupa and Rocha-Pinto2017).

We only proceeded with a simulation when the resulting tidal radius was smaller than its Jacobi radius at the cluster’s initial galactic radius. For computational efficiency, the Jacobi radius $r_J$ was calculated assuming a circular galactocentric orbit and a flat rotation curve, giving

(17) \begin{equation} r_J = \left(\frac{{\rm GM}_cR_g^2}{2v_c^2}\right)^{1/3},\end{equation}

where $M_c$ is the cluster mass, $R_g$ is the galactic centred radius and $v_c$ is the circular velocity of the Galaxy at that radius.

This condition ensures that the cluster initially lies within its tidal boundary in the Galactic potential, preventing stars at large radii from being immediately stripped by the galactic tidal field. Without this, the input parameters would not represent a physically stable initial cluster configuration, affecting the relationships between input and output parameters that the machine learning models are trained on. We also assume that all neutron stars and black holes are ejected from the cluster due to the large velocity kicks at their creation and the relatively low escape velocity of the investigated clusters (Pavlík et al. Reference Pavlík, Jeřábková, Kroupa and Baumgardt2018).

Initial conditions were sampled from a 7-dimensional Latin hypercube spanning a broad range chosen to encompass all physically reasonable initial conditions, ensuring adequate coverage of parameter space for emulator training. If the best fitting parameters were later found to be in a region with a low density of training data, we ran additional simulations in that region of parameter space. Once all simulations were completed, we applied the same $\chi^2$ filtering described in Equations 1 and 2 on the final state of the simulated clusters so that the present day values in the simulations could be compared with observations.

We then measured the final number of systems, binary fraction and 2D half mass radius of each simulated cluster as well as the mass function using the same break points as the observed clusters. The surface density of each cluster was also measured by binning stars into 10 radial bins with an equal number of stars in each, and finding

(18) \begin{equation} \rho_{{\rm surf},i} = \frac{N_{\rm star},i}{\pi \left(r_i^2-r_{i-1}^2 \right)}\end{equation}

for each bin i where $r_i$ is the outer radius of bin i.

We then found the surface density of each simulation’s final state using the same radial bins as the observed clusters. Thus, in total we have 16 test statistics that we can use to compare the simulated clusters to the observations. The aim of this analysis is to create an emulator that will take any combination of the 7 simulation input parameters and give a prediction of these 16 test statistics without needing to process an entire N-body simulation of the cluster.

Table 5.

Best fitting initial cluster parameters. Reported values include the inferred initial number of stellar systems ( $N_{\text{sys,ini}}$ ), initial binary fraction ( $f_\mathrm{bin,ini}$ ), initial 3D half-mass radius ( $r_{h,{\rm ini}}$ ), the system IMF slopes for the low, intermediate and high mass ranges ( $\alpha_{l,{\rm sys}}$ , $\alpha_{m,{\rm sys}}$ , $\alpha_{h,{\rm sys}}$ ) and the stellar IMF slopes for the low, intermediate and high mass ranges ( $\alpha_{l,{\rm star}}$ , $\alpha_{m,{\rm star}}$ , $\alpha_{h,{\rm star}}$ ). Uncertainties correspond to 1 $\sigma$ errors.

Table 5. Long description

A table comparing initial cluster parameters for different star clusters. The table has 4 columns and 4 rows, including the header row. The columns are labeled as Cluster, N sys ini, f bin ini, r h ini (pc), alpha 1 sys, alpha m sys, alpha h sys, alpha 1 star, alpha m star, and alpha h star. The rows are labeled with the names of the star clusters: Alpha Persei, Pleiades, and Praesepe. Each row provides specific values for the initial number of stellar systems, initial binary fraction, initial 3D half-mass radius, and the system and stellar IMF slopes for the low, intermediate, and high mass ranges. The values include uncertainties corresponding to 1 sigma errors. For example, Alpha Persei has an initial number of stellar systems of 1127 plus or minus 30, an initial binary fraction of 0.346 plus or minus 0.019, and an initial 3D half-mass radius of 5.63 plus or minus 0.35 parsecs. The table provides a detailed comparison of these parameters across the different star clusters.

Figure 11.

The derived highest likelihood stellar IMF (red) for each cluster compared to its PDMF (black). The mass function $\xi(\log m)$ is scaled such that the curves reflect simulated stellar losses of $\approx15\%$ , 9% and 49% for Alpha Persei, Pleiades, and Praesepe, respectively. We see that Praesepe is more evolved and thus has a larger number of low mass stars that have been stripped over time, whereas Alpha Persei and Pleiades are younger and have not gone through as much significant dynamical evolution.

To create this emulator and find the optimal method of emulation, we use the python package AUTOEMULATE (Stoffel et al. Reference Stoffel2025). Using this tool, we compare various models to find the method which maximises the $R^2$ statistic when applied to a test subset of data. From this analysis, we found that using a correlated Gaussian process was optimal for modelling the relationships between the input parameters and test statistics.

On initial inspection, we found that the concentration parameter c had little significant correlation with any observable test statistic and was best tuned so that the initial cluster had a tidal radius approximately equal to the Jacobi radius for a given initial half mass radius. Thus, we omitted the input variable c in the emulator training model to reduce complexity in the emulator and the subsequent solver.

5.2 Parameter fitting

To find the initial cluster conditions that are most likely to recreate observations, we used the inbuilt Bayesian modelling in AUTOEMULATE. More specifically, we assume that all observables have Gaussian errors and implement a Hamiltonian Monte Carlo (HMC) No-U-turn sampler (NUTS) algorithm to sample the posterior distribution. We then report the median values of the sample as the best fitting values with uncertainties given by the 16th and 84th percentiles. The resulting posterior distributions can be seen in Figure B1 and the best fitting values in Table 5.

Due to the way the cluster is created for the simulations, we solved for the maximum likelihood estimators of the initial system IMF and binary fraction rather than the initial stellar IMF. Similarly to Section 4.2, we create a large synthetic cluster to estimate the corresponding stellar IMF of each cluster. Using this method, we find the true IMF of the three clusters in Table 5 and Figure 11.

We see that Alpha Persei and Pleiades undergo minimal dynamical evolution which is consistent with expectations due to their relatively low age with respect to their relaxation times. This is in contrast to Praesepe, which is seen to have lost a significant fraction of lower mass stars below $\approx0.6\,\mathrm{M}_\odot$ due to mass segregation and prolonged tidal stripping. Interestingly, we find that the high mass IMF of each cluster is steeper than the observed PDMF by about 0.2. This is likely driven more by the dependence of the final cluster state’s half mass radius on the initial high mass function slope as a larger number of high mass stars will cause additional dynamical heating as well as drive stronger expansion early in the cluster’s evolution from the shallower gravitational potential caused by early, high mass star evolution. Thus, the slope of the IMF at high masses is constrained on both the measurement of the high mass PDMF and the present day half mass radius.

5.3 Validation of initial conditions

We verified the best fitting parameters of each cluster using N-body simulations and comparing the final states to the observed clusters. Since every created cluster and simulation will have random variations, especially with the smaller sample of stars in open clusters, we ran at least 20 simulations with the best fitting parameters for each cluster and chose the run with the best agreement with observations.

We find that the simulations are able to reproduce the observed number of stars, binary fraction, half mass radius and three mass function slopes within $2.5\sigma$ , with most parameters agreeing within $1\sigma$ . In addition, while we do not include any indicators of mass segregation into the optimised likelihood function, we find good agreement in the degree of mass segregation between observations and simulations. A Kolmogorov–Smirnov (KS) test between each of the observed and simulated radial distributions found no statistically significant differences, with p-values ranging between 0.15 and 0.77. As shown in Figure 13, the degree of mass segregation for each cluster is related to their age with respect to their relaxation times. The difference in the observed low mass and high mass radial distributions in Alpha Persei shows no statistically significant difference, with $p=0.30$ , whereas Pleiades and Praesepe show significant differences in their stellar radial distributions with $p\lt 10^{-5}$ and $p\lt 10^{-15}$ respectively.

Figure 12.

The observed surface density profile of Praesepe (black) and the surface density profile of a simulation using the best fitting parameters (red dashed). There is good agreement between the two curves, showing the validity of the derived parameters. We also show the best fitting Plummer model for the observed surface density (blue dot-dashed) with Plummer radius $a=4.26\pm0.11$ and $N_\mathrm{star}=1367\pm40$ .

Figure 13.

Cumulative distribution of stars in each cluster with respect to radius (black). Cluster members are split into two groups, showing the different radial distribution of high mass stars $\gt 1\,\text{M}_\odot$ and low mass stars $\lt 1\,\text{M}_\odot$ which indicates the presence of mass segregation. We also show the radial distribution of simulations using the best fitting parameters (red) which are in good agreement, despite mass segregation not being used in the goodness of fit estimation.

We also find good agreement with the simulated and observed surface density profiles for each cluster, with Praesepe shown in Figure 12. Röser & Schilbach (Reference Röser and Schilbach2019) found that the surface density of Praesepe is well fitted with a Plummer model (Plummer Reference Plummer1911) with core radius $r_\mathrm{co}=3.7\,\mathrm{pc}$ . We fit a 2D Plummer model to the surface number density of each cluster, given by

(19) \begin{equation} \rho_\mathrm{surf}(r) = \frac{N_\mathrm{star}}{\pi a^2}\left(1+\frac{r^2}{a^2}\right)^{-2},\end{equation}

where $N_\mathrm{star}$ is the total number of stars in the cluster and a is the Plummer radius.

We fit both $N_\mathrm{star}$ and a as free parameters for each cluster using a non-linear least squares fit. For Alpha Persei, we find the best fitting Plummer model is given by $a=7.17\pm0.20$ pc and $N_\mathrm{star}=1\,105\pm35$ which has a reduced $\chi^2$ value of $\chi^2_\mathrm{red}=2.35$ . The surface density profile of Pleiades is fitted by $a=4.13\pm0.10$ pc and $N_\mathrm{star}=1\,646\pm42$ with $\chi^2_\mathrm{red}=8.48$ and Praesepe is described best by a model with $a=4.26\pm0.11$ pc and $N_\mathrm{star}=1\,367\pm40$ with $\chi^2_\mathrm{red}=3.57$ . While all three models give $\chi^2_\mathrm{red}\gt 1$ , the Plummer model still provides a reasonable approximation of the surface density profile of each cluster using a simple analytical model. This can be useful for obtaining a first-order estimate of the structural scale of each cluster and for providing a simple analytic description of the observed stellar distributions.

The agreement between these forward simulation runs and cluster observations help demonstrate that these solutions are valid initial conditions of the three investigated clusters. These values, however, are not necessarily the only possible initial state of the cluster. Wang & Jerabkova (Reference Wang and Jerabkova2021) demonstrate that small changes in the number of high mass OB stars have significant effects on the evolution of a cluster over the span of a few 100 Myrs. If we assume that the IMF is randomly sampled during star formation and acts as a probability distribution function, similar to the sampling methods of this work, then stochastic variations in the number of high-mass stars can significantly alter the subsequent dynamical evolution and present-day state of a cluster compared to one drawn from the same IMF but containing fewer massive members. There are also additional variations introduced due to random dynamical encounters as well as the random initialisation of binaries in each N-body simulation. We expect all of these stochastic features to introduce scatter in the training data which the emulator models are trained on, however, the degree of variations should be reflected by uncertainties of the emulators predicted values. These uncertainties will then be reflected by a wider spread of data when sampling the posterior distribution and thus a larger uncertainty in the predicted initial parameters.

Thus, in this work we have found gas-free initial conditions of the present-day cluster state that have the highest likelihood of reproducing present-day observations, but will have variations with each realisation of the cluster. While the true initial conditions may vary, we show that this method is able to accurately recreate observations and provides a useful initial state for simulation work that wishes to study the evolution of a cluster under different conditions.

5.4 Comparison to other works

While a direct comparison between different IMF measurements is difficult due to the different mass ranges which the mass function slopes are calculated over, we compare our results to other IMF determinations of the same three clusters as well as other system types for similar mass ranges. As shown in Figure 14, our IMF values agree with most other determinations of cluster IMFs over similar mass ranges with only a few exceptions. The largest discrepancies are likely due to the significantly different mass ranges which the power law is calculated over. In many studies of star cluster IMFs, the adopted break points for the multi-stage power law are ones that give a global best fit, rather than being specific to each cluster. While this fits into the framework of a universal cluster IMF, it can result in a IMF that does not well describe all of the features in a particular cluster and can smooth out unique turnover behaviour, resulting in what appears to be a universal average mass function slope that is not necessarily a very good fit for a number of clusters.

Figure 14.

Comparison of IMF measurements of Alpha Persei, Pleiades, Praesepe, and other stellar bodies. The mass range in M $_\odot$ which the value was calculated over is noted underneath each data point. The red vertical line and shaded region depict the mean IMF slope values for our three measured clusters with uncertainty. The grey vertical line corresponds to the slope values from the canonical Kroupa (Reference Kroupa2001) IMF.

In our approach, we determine the break points for each cluster independently and find that a high mass function break of $m_{x2}=1\,\mathrm{M}_\odot$ is within error bars for all three clusters, consistent with the majority of other IMF determination studies. The low-intermediate mass break point $m_{x1}$ , however, is less defined and has greater variation between the three clusters as well as with other studies of the same clusters. While there could be a real variation of turn over behaviour at low masses between clusters, the specific values of $m_{x1}$ we find are highly dominated by the choice of mass-luminosity relation or isochrone, meaning we cannot confidently determine any physical relationships.

We find nearly all our individual cluster measurements to be in general agreement with other determinations over similar mass ranges. For Alpha Persei, Sheikhi et al. (Reference Sheikhi2016) find that $\alpha_{l}=0.89\pm0.11$ and $\alpha_{h}=2.37\pm0.09$ using the PPMXL and Wide-field Infrared Survey Explorer catalogues, correcting for unresolved binaries, and Lodieu et al. (Reference Lodieu, Deacon, Hambly and Boudreault2012) find that $\alpha_{l}=0.64\pm0.08$ using UKIDSS data. These values are in good agreement with our findings of $\alpha_{l}=0.98\pm0.23$ , $\alpha_{m}=1.64\pm0.25$ and $\alpha_{h}=2.60\pm0.19$ . The slight difference in the high-mass function slope likely arises from our use of two power laws over 0.30–0.95 and 0.95–5.55 M $_\odot$ , rather than a single power law over 0.62–4.68 M $_\odot$ as adopted by Sheikhi et al. (Reference Sheikhi2016). Our IMF slopes for Pleiades in the intermediate and high mass ranges also agree with works such as Niu et al. (Reference Niu, Wang and Fu2020) and Sollima (Reference Sollima2019) who find $\alpha_{m}=1.97\pm0.22,1.6\pm0.2$ and $\alpha_{h}=3.09\pm0.18,3.4\pm0.3$ respectively using LAMOST and Gaia DR2 observations.

Lastly, we also find good agreement with previous determinations for Praesepe for the high and intermediate MF slopes. Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013) found that $\alpha_{h}=2.88\pm0.22$ which agrees well with our measured value of $\alpha_{h}=2.86\pm0.36$ . Our measured intermediate IMF slope of $\alpha_{m}=2.23\pm0.24$ also agrees well with the findings of Ebrahimi et al. (Reference Ebrahimi, Sollima and Haghi2022) who determined $\alpha_{m}=2.16\pm0.13$ , although these values are all calculated over slightly different mass ranges which may lead to slightly different values. For the low mass IMF slope, however, we find a much steeper value compared to published values. Khalaj & Baumgardt (Reference Khalaj and Baumgardt2013) found that $\alpha_{l}=0.85\pm0.10$ , and Baker et al. (Reference Baker2010) determined that $\alpha_{l}=1.10\pm0.37$ using near infrared survey data, whereas we find a higher value of $\alpha_{l}=1.93\pm0.24$ . The main contributing factor to this discrepancy is likely our inclusion of dynamical evolution which would include many additional low mass stars compared to the present-day state of the cluster. Moreover, there is additional error in the exact slope value due to the uncertainties of the mass conversion models for these low masses as previously discussed.

In addition, we also find that simulations of Praesepe which use the best fitting parameters can qualitatively reproduce the tidal tails observed by Röser & Schilbach (Reference Röser and Schilbach2019). In our simulations, we identify a total of 700 stars in the tidal tails outside of two tidal radii and within a radius of 165 pc, matching the observational bounds of Röser & Schilbach (Reference Röser and Schilbach2019) who used Gaia DR2 data. Of these stars, we find 339 in the leading tail and 361 in the trailing tail. This is broadly consistent with the 389 tail stars identified by Röser & Schilbach (Reference Röser and Schilbach2019), particularly considering that observations of the trailing arm are likely incomplete for fainter stars due to its larger distance from the Sun. The overall agreement suggests that our simulations are able to capture the dominant dynamical processes responsible for the formation of these extended tidal features.

We note, however, that our simulations assume gas-free initial conditions. As shown by Weis et al. (Reference Weis2025), the inclusion of instantaneous gas expulsion at cluster birth can remove up to $\approx75\%$ of the initial cluster mass, implying a significantly larger initial stellar population than considered here. In the work of Weis et al. (Reference Weis2025), this leads to a considerably different inferred IMF of Praesepe, with a very steep high mass slope ( $\alpha_{h}=4.1$ ) and much fewer low mass stars ( $\alpha_{l}=0.4$ ). However, their analysis does not consider unresolved binaries which, as discussed in Section 4, can bias the inferred mass function. As such, the differences between these findings reflect both differing initial assumptions and modelling choices. We also note that a range of similar dynamical modelling approaches have previously been applied to star clusters. For example, Kroupa et al. (Reference Kroupa, Aarseth and Hurley2001) modelled the early gas-embedded evolution of a Pleiades-like cluster from Orion Nebula Cluster-like initial conditions, while Haghi et al. (Reference Haghi, Zonoozi, Kroupa, Banerjee and Baumgardt2015) show the effect of gas expulsion on the MF in the first 100 Myr of the evolution of globular clusters. However, Haghi et al. (Reference Haghi, Zonoozi, Kroupa, Banerjee and Baumgardt2015) found that gas expulsion alone does not strongly alter the inferred IMF unless strong primordial mass segregation is also present, which we do not include in our models.

To compare our results to the mass functions of other stellar groups, we found the mean IMF slope of the three open clusters for each mass range. We also determined the standard deviation of each IMF slope similar to the formalism in Dib et al. (Reference Dib, Schmeja and Hony2017), where the IMF has cluster-to-cluster variations that can be described by Gaussian distributions. From the three clusters investigated, we find the average IMF to be given by $\bar\alpha_l=1.24\pm0.29$ , $\bar\alpha_m=1.72\pm0.09$ and $\bar\alpha_{h}=2.98\pm0.22$ with dispersions of $\sigma_{l}=0.43\pm0.24$ , $\sigma_{m}=0.00\pm0.14$ , $\sigma_{h}=0.29\pm0.16$ and break masses of $m_{x1}=0.25\pm0.02\,\mathrm{M}_\odot$ and $m_{x2}=0.96\pm0.14\,\mathrm{M}_\odot$ .

This average open cluster IMF is shown in Figure 15 against the canonical field IMF from Kroupa (Reference Kroupa2001) and the measured average IMF of globular clusters as determined in Baumgardt et al. (Reference Baumgardt, Hénault-Brunet, Dickson and Sollima2023), where all functions are normalised so that the area under the IMF is unity. In comparison to both functions, the open cluster IMF depicts a smaller population of high mass stars $\gt 1\,\mathrm{M}_\odot$ . While the values of the break points and slopes are different to those quoted in Kroupa (Reference Kroupa2001), the turnover shape of the IMF for low and intermediate mass stars below $\approx1\,\mathrm{M}_\odot$ is qualitatively similar. This is in contrast to the bottom-light IMF form observed in globular clusters which have a significantly smaller proportion of low mass stars $\lt 0.3\,\mathrm{M}_\odot$ and larger proportion of intermediate and high mass stars. The apparent deficit high mass stars in the investigated metal-rich open clusters compared to metal-poor, denser globular clusters is in qualitative agreement with previous studies that suggest denser star forming environments produce increasingly top-heavy IMFs (Marks et al. Reference Marks, Kroupa, Dabringhausen and Pawlowski2012; Dabringhausen et al. Reference Dabringhausen, Kroupa, Pflamm-Altenburg and Mieske2012).

Figure 15.

Average open cluster initial mass function (IMF) (black solid line) compared to the canonical IMF from Kroupa (Reference Kroupa2001) (red dashed) and the average globular cluster IMF from Baumgardt et al. (Reference Baumgardt, Hénault-Brunet, Dickson and Sollima2023) (blue dot-dashed). All three curves are normalised such that the area under the curve is unity. The open cluster IMF closely follows the canonical shape but exhibits a flattening at intermediate masses and a steepening at high masses, implying a relative deficiency of massive stars. Relative to the globular cluster IMF, it contains a larger fraction of low-mass stars and fewer high-mass stars.