The Present-day Mass Function of Star Clusters in the Solar Neighborhood

This work analyzes the present-day mass function (PDMF) of 93 star clusters utilizing Gaia Data Release 3 data, with membership determined by the StarGo machine-learning algorithm. The impact of unresolved binary systems on mass estimation is rigorously assessed, adopting three mass ratio profiles for correction. The PDMF is characterized by the power-law index, α, derived through a robust maximum likelihood method that avoids biases associated with data binning. The value of α for stars between the completeness limited mass of Gaia (with a mean 0.3 M ⊙ for our cluster samples) and 2 M ⊙ exhibits stability for clusters younger than 200 Myr, decreasing for older clusters, particularly when considering stars within the half-mass radius. The PDMF of these star clusters is consistent with a dynamically evolved Kroupa initial mass function via the loss of low-mass stars. Cluster morphology shows a correlation with α, as α values exhibit a decreasing trend from filamentary to tidal-tail clusters, mirroring the sequence of increasing cluster age. The dependence of α on the total cluster mass is weak, with a subtle increase for higher-mass clusters, especially outside the half-mass radius. We do not observe a correlation between α and the mean metallicity of the clusters. Younger clusters have lower metallicity compared to their older counterparts, which indicates that the older clusters might have migrated to the solar neighborhood from the inner disk. A comparison with numerical models incorporating a black hole population suggests the need for observations of distant, older, massive open clusters to determine whether or not they contain black holes.


INTRODUCTION
The concept of initial mass function (IMF), describing the distribution of stellar masses at birth, was first proposed by Salpeter (1955).The IMF is mathematically formulated as a power-law function, where A is a normalization coefficient, m is the stellar mass, and ξ(m) is a function that is either the number of stars or the normalized fraction of stars within the mass range m + dm.α is the power-law index of the Xiaoying.Pang@xjtlu.edu.cnmass function, and α = 2.35 is referred to as "Salpeter slope".
The IMF holds significant importance in comprehending numerous astronomical phenomena, including the formation of the first stars (Bromm et al. 2009), the evolution and formation of galaxies (Calura & Menci 2009;Yan et al. 2021), and the precise determination of the absolute star formation rate (Jeřábková et al. 2018a;Aoyama et al. 2023).Furthermore, the IMF plays a vital role in the theoretical framework concerning star formation, representing the ultimate outcome of molecular cloud contraction and fragmentation processes (Krumholz 2014).
Accumulating evidence suggests a deviation from the "Salpeter slope" in the mass function of lower-mass stars Pang et al. in the Galactic disc (Kroupa et al. 1993;Kroupa 2001;Chabrier 2003).In addressing this, Kroupa et al. (1993) proposed an alternative IMF, characterized by a multisegment power law.This "Kroupa IMF" exhibits a shallower power-law index at lower masses (in contrast to the Salpeter value, which is obtained from higher mass stars).
Numerous endeavors have been undertaken to determine the IMF in the Milky Way.However, the Galactic field star population represents a diverse assembly of stars characterized by different ages and metallicities.Prior examinations of the field IMF frequently limit their focus to samples of low-mass stars to mitigate evolutionary effects over the age of the universe (Kroupa et al. 1993;Li et al. 2023).An accurate star formation history estimation of our galaxy, which degenerates with the estimation of the IMF for intermediate-mass and massive stars, remains challenging (Mor et al. 2019).
Young star clusters are optimal targets for determining the IMF.They contain a statistically significant number of stars and share similar characteristics, such as distance, age, and chemical composition.Moreover, a significant portion of stars in the Galactic field have their origin in unbound stellar groups that quickly dissipate after gas expulsion (Lada & Lada 2003).The filamentary-type stellar group identified in Pang et al. (2022a) is one of these kinds, which is formed in environments with low star formation efficiency (Kruijssen 2012).Analyzing the mass function of such groups holds the potential to enhance our understanding of the role of dissolution in the evolution of mass distribution in these young stellar groups.
A major limitation in studying the IMF of star clusters arises from their dynamical evolution over time.This evolution leads to a systematic depletion of low-mass stars through evaporation, and the initial loss of certain massive stars induced by dynamical ejection resulting from binary encounters, and subsequently due to stellar evolution.This process results in a stellar mass function that changes with cluster age (Portegies Zwart et al. 2010;Vesperini & Heggie 1997;Baumgardt & Makino 2003).Therefore, accounting for the dynamical evolution of star clusters becomes essential when attempting to infer the IMF of a cluster from its present-day mass function (PDMF).Nearly all open clusters exhibit mass functions that resemble, or are consistent with, a dynamically evolved Kroupa/Chabrier-type IMF (Bastian et al. 2010).
When correcting the PDMF of star clusters for evolution, to obtain the IMF, all the aforementioned factors introduce uncertainties.Therefore, many measurements are typically presented in the form of the PDMF, with authors not attempting to deduce the IMF to avoid introducing additional uncertainties.
To obtain an accurate description of the PDMF, star clusters with reliable memberships and high-quality photometry observations are essential.A robust membership avoids significant field contamination, preserving the anticipated PDMF.Precise photometry enables the accurate derivation of individual stellar masses.There are several additional factors that we need to consider when deriving the PDMF.The presence of unresolved binary stars should be accounted for, to avoid the introduction of a slight shift in the observed PDMF peak toward higher masses (Lee et al. 2020).Simultaneously, the observational completeness of low-mass stars affects PDMF.For example, the 21 mag limit in the G band of Gaia (Gaia Collaboration et al. 2022) introduces limitations in detecting the lowest mass stars in star clusters.In particular, discerning evaporated low-mass stars situated in the periphery or in the tidal tail of a cluster presents a considerable challenge.The mass segregation effect biases the measured stellar population in the inner radius to massive stars (Gouliermis et al. 2004;Allison et al. 2009;Pang et al. 2013).Despite utilizing spacebased observatories like Gaia, the pronounced crowding of stars in the central region of distant clusters remains a formidable challenge for the observation of faint stars (Baumgardt et al. 2023).
To satisfy the aforementioned criteria, and avoid the crowding effect of distant clusters, this work we focuses on star clusters about 500 pc of the solar neighborhood.We are motivated to explore variations in the PDMF among star clusters of different ages and environments, aiming to reveal the PDMF of various star cluster types, to examine the influence of dynamical evolution on the mass distribution, and to identify the key parameters that affect the mass distribution in star clusters.
This paper is organized as follows.In Section 2, we introduce the dataset and membership of star clusters.The PDMF is derived in Section 3. In Section 3.1, we adopt three different mass ratio profiles and correct the observed stellar mass for binaries.The power-law index of the PDMF is obtained via the maximum likelihood method in Section 3.2.In Section 3.3, we discuss the evolution of maximum stellar mass in star clusters.The dependence of the PDMF on cluster parameters is investigated in Section 4. In Section 4.1, we present how dynamical evolution affects the PDMF.We then investigate the relation between PDMF and the cluster chemical abundance in Section 4.2.In Section 4.3, we compare the PDMF to numerical models including black holes.Finally, we provide a brief summary of our findings in Section 5.

MEMBERSHIP FROM GAIA DATA
In this study, we investigate a total of 93 star clusters.The identification of member stars of 92 of the target star clusters was carried out in Pang et al. (2021a,b, 2022b,c), and Li et al. (2021).These studies established stringent criteria, requiring member stars to have parallax and photometric measurements with uncertainties within a 10% threshold.The selection of star cluster members is executed using the machine learning algorithm StarGo (Yuan et al. 2018) based on Gaia EDR 3 and DR 3 data (Gaia Collaboration et al. 2021).Member stars are chosen with a contamination rate of at most 5%, resulting in a corresponding membership probability of 95%.The precision of the photometry has been improved for sources with magnitudes fainter than G = 13 mag following the correction of G-band photometry in Gaia DR 3 by Riello et al. (2021).To enhance the accuracy of photometric and kinematic data for each star, we cross-match the members with Gaia DR 3 (Gaia Collaboration et al. 2022), and our analysis relies on the data from Gaia DR 3. We identify members of one old open cluster NGC 6991 (> 1 Gyr, Cantat-Gaudin et al. 2020) in this work.The membership of this cluster was determined with the same method described above, using Gaia DR 3 data (detailed in Appendix A).In total, we have 93 star clusters in our sample for PDMF analysis.The basic parameters for these 93 clusters are presented in Table 1.  ). αPDMF is the power-law index of the PDMF for stars with a mass between m lower and 2 M⊙.M cl is the observed total stellar mass above the completeness mass limit of a star cluster.mmax is the observed maximum stellar mass in each cluster.m lower is the completeness limited mass of each cluster from Gaia DR 3 data.Columns 7 and 8 are mean metallicities of target clusters derived from Gaia DR 3 XP spectra (Li et al. 2023b) and APOGEE, respectively.The error is the dispersion of the metallicity of each cluster.The last column is the morphological type taken from Pang et al. (2022b).

Binary correction
To construct the PDMF of the target clusters, it is necessary to first compute the masses of individual member stars.Previous studies (Pang et al. 2021a(Pang et al. ,b, 2022b,c;,c;Li et al. 2021) have obtained the best-fitting isochrones for the 92 star clusters by assessing the locations of member stars in the color-magnitude diagram (CMD).The same procedure is carried out for NGC 6991 in this work (Figure A.1 in Appendix A).Previous studies directly used the isochrone provided by PARSEC.This time we enrich the data points in each isochrone by a factor of four through interpolation.Following previous studies (Liu & Pang 2019), we search for the nearest neighbour of each observed star on the isochrone using the KD-tree method (McMillan et al. 2007), and assign the stellar mass of the nearest neighbor to the observed star as its stellar mass, m 0 .
The uncertainty in the stellar mass m 0 is primarily determined by the photometric uncertainty of Gaia DR3 (Gaia Collaboration et al. 2022), which may result in variations in the locations of cluster members in the CMD, and variations in the estimation of the binary properties (Pang et al. 2023).To determine the impact of these photometric uncertainties, we generate simulated observations.We assign each member a new magnitude drawn from a Gaussian distribution centered at the observed value and with a dispersion equal to the photometric uncertainty of the individual star.With the updated location of each member in the CMD, we re-peat the aforementioned procedure, and assign the stellar mass, m ′ 0 , that is equal the nearest point on the best-fitting isochrone of each cluster.This procedure is repeated 1000 times.The difference between the m 0 and m ′ 0 is the adopted uncertainty of the stellar mass in this work.The average uncertainty in the stellar mass is 0.01 M ⊙ .
The stellar mass, m 0 , obtained in the manner described above is determined under the assumption that observed stars are single stars.However, as shown by Pang et al. (2023), the unresolved binary fraction (for mass ratio q >0.4) among 85 star clusters ranges from 5% up to 40% (Figure 2 (a) in Pang et al. 2023).Note that these values are lower limits of the total unresolved binary fraction.Consequently, in order to reliably quantify the PDMF of the clusters in our sample, we must correct individual stellar masses for binary effects.
If an observed star is indeed an unresolved binary system, its total mass is under-estimated when using m 0 .Our objective is to identify the corresponding primary and secondary masses for potential binary systems, which will be incorporated into the PDMF. Lee et al. (2020) demonstrated that the influence of binarity on the mass function is sensitive to both the binary fraction and the mass ratio of the systems involved.
To achieve this goal, we implement three mass ratio profiles in the binary correction, namely the uniform distribution, the profile of Kouwenhoven et al. (2007), and the profile of Fisher et al. (2005)  ter.Each point on the isochrone is treated as a single star, with its corresponding mass m ′ designated as the primary mass.A mass ratio q is randomly generated from one of the three the mass ratio profiles.A companion star of mass qm ′ , is assigned from the isochrone.This process is repeated 100 times.For the low-mass primaries approaching the limit of 0.08 M ⊙ , their secondary mass qm ′ is not available in the isochrone when q is small.Consequently, the adopted q values for these low-mass stars are preferentially higher than what is expected from the three adopted mass ratio distributions (see Kouwenhoven et al. 2010, for details).We combine the luminosity of both primary and secondary stars in the BP , RP , and G bands to compute the magnitude for a synthetic unresolved binary system in these filter bands.Three populations of synthetic unresolved binary systems are generated following the aforementioned three mass ratio distributions.Given that the mean completeness mass limit of Gaia data for our cluster samples is 0.3 M ⊙ , exceeding 0.08 M ⊙ , the uncertainty associated with low-mass synthetic unresolved binaries featuring a primary component of m ′ ∼ 0.08 M ⊙ , does not impact on our findings (see Section 3.2).
We assign a binarity probability p for each observed star based on the empirical relation between binarity frequency and primary mass m ′ (the crosses in Figure 2) taken from Offner et al. (2022).When adopting this relation, we assume that the correlation in question remains unaltered across varying cluster ages and the potential influence of dynamical evolution is neglected.We note that a higher initial binary fraction for low-mass stars has been suggested which decreases over time to the present-day binary fraction due to dynamical interactions (Kroupa 1995a,b;Thies et al. 2015).However, this dynamical evolution primarily results in the disruption of wide binary systems, while the impact on the hard binary systems, which are unresolved in our observations, is negligible.Given the substantial discrepancies in the observed and simulated dependencies of p on m ′ , we consider the adoption from Offner et al. (2022) as an adequate choice.By assuming the stellar mass m 0 of the observed stars as the primary mass, we interpolate the value of p for an observed star based on this relation (orange curve in Figure 2).Utilizing this probability, each observed star is assigned either as an unresolved binary system or a single star.In the case of an unresolved binary system, we identify the synthetic unresolved binary system from the binary population, whose position is closest to the observed star in the CMD.The corresponding primary and secondary masses of the binary system are considered in computing the PDMF (instead of m 0 ).This process is repeated 1000 times, during which the observed stars are treated as unresolved binary systems in p × 1000 instances, and as single stars in (1 − p) × 1000 instances.Each process involves the determination of the PDMF.

Power-law index α of the PDMF
The Kroupa IMF (Kroupa 2001) has a power-law index equal to 2.3 for m > 0.5 M ⊙ (see section 3.2 in Yan et al. 2023a), when the influence of stellar density and metallicity on the IMF for initial cluster conditions is neglected (Jeřábková et al. 2018b).When determining the shape of the PDMF over a small stellar mass range, it is in many cases practical to adopt a single power law.This approach avoids the complications of the multiple parameters (the break point, and the two power indices) of a PDMF with two power law segments.This approach is commonly applied in star counting studies with a stellar mass detection limit similar to our completeness limit 0.3 M ⊙ (Geha et al. 2013;Gennaro et al. 2018;Yang et al. 2024).Therefore, we model the PDMF for our target clusters as a single power-law profile for the masses above the typical completeness limit 0.3 M ⊙ in our sample.
To determine the PDMF, we need to fit two parameters in Equation 1: the normalization parameter A, and the power-law index α.These parameters are sensitive to the approach how the data are binned (Maíz Apellániz & Úbeda 2005).To mitigate the bias arising from data binning, we adopt the maximum likelihood method to derive these coefficients.
To establish the likelihood function for a group of stellar masses, we assume that each mass is independently measured.Additionally, we postulate that the probability of obtaining a specific mass measurement is the mass function's value at that particular mass in Equation 1.
A sample of stellar masses of each cluster (m i , where i ranges from 1 to N , and N represents the total number of members in each cluster) is taken after the binary correction.The likelihood function for this sample can be formulated as the product of the individual probabilities associated with each m i : (2) Therefore, the log-likelihood of Equation 2 can be written as Since N is the integral of the PDMF from the lower stellar mass limit (m lower ) to the upper stellar mass limit(m upper ), Here we express the parameter A as a function of total member number N , m lower , m upper and power-law index α in Equation ( 4).Given the incompleteness of the Gaia data, our member stars are subject to incompleteness beyond a specific magnitude.To avoid the bias introduced by incompleteness, we set the m lower for each cluster as the completeness limited mass (dashed vertical red lines in Figure 3), whose value is presented in column 6 in Table 1.The maximum stellar mass in each cluster varies significantly.Old clusters lack high-mass stars as a result of stellar evolution.To ensure a consistent mass range in the analysis of all target clusters, we only estimate the PDMF slopes for low-mass stars and adopt m upper = 2 M ⊙ (dotted vertical red lines in Figure 3) for all clusters.
In this work, we utilized the NumPyro package1 , for performing Markov Chain Monte Carlo (MCMC) sampling.Flat priors are assigned to the parameters A and α, and a total of 1000 warm-up iterations are performed for each cluster, followed by 1000 sampling iterations.After this sampling process, the mean value of α is considered as the most probable power-law index of the PDMF between the completeness-limited mass to 2 M ⊙ for each cluster.The mean standard deviation of α of all 1000 repeated processes is treated as its uncertainty.The value of the α and its uncertainty for each cluster is presented in column 3 of Table 1.
In Figure 3 we present two examples of the PDMFs: Pleiades (left panels) and Praesepe (right panels).The most probable PDMFs (between completeness limit and 2 M ⊙ ) are indicated with orange dashed lines.The histograms for the values of α values for all 93 target clusters are shown in Figure 4.
The peak at α ≈ 2.0, agrees with the value of the "Kroupa IMF" (Kroupa 2001), computed for the mass range of 0.3 M ⊙ to 2 M ⊙ : α = 2.04.The distribution of α is nearly identical for the three mass ratio distributions adopted in the binary correction, indicating a negligible difference induced by the choice of the mass ratio profiles.Consequently, we only use the α value obtained from the uniform mass ratio distribution for subsequent analysis in the following sections.

Maximum stellar mass in the cluster
In Figure 5 (a) we present the dependence of the observed maximum stellar mass, m max (column 5 in Table 1), on the observed total stellar mass above the completeness mass limit, M cl (column 4 in Table 1).Each data point represents one cluster, with its color indicat- ing the cluster age.While the mass estimation method introduced in the previous sections is suitable for a stellar population, the uncertainties arising from dust extinction, binary probability, binary evolution, merger probability, and dynamical ejection become inherently irreducible through statistical arguments when assessing the mass of a single star.Addressing these factors typically yields a mass uncertainty exceeding 0.3 dex (Yan et al. 2023b).Conducting a detailed mass uncertainty analysis through dynamical simulations of star clusters is complicated and lies beyond the scope of the current study.Consequently, we adopt the largest m 0 as m max without accompanying uncertainty estimation.The findings in this section should be interpreted as suggestive rather than conclusive.
As can be seen from Figure 5 (a), older clusters (orange to yellow dots) have a smaller maximum observed stellar mass, which can be primarily attributed to stellar evolution.This trend is more prominent in panel (b), which shows the relation between m max and the cluster age.A positive correlation is observed between m max and M cl (grey curve, panel (a)), indicating that more massive young clusters generally host stars with a higher m max .This finding aligns with the observed m ini max -M ecl relation reported in previous studies (Weidner et al. 2010;Yan et al. 2017;Yan et al. 2023b), in which m ini max is the initial maximum stellar mass in a young star cluster, and M ecl is the embedded cluster mass.To further illustrate this, we compare the observed m max values with the theoretical stellar maximum mass limit ( thin dashed lines).Adopting the stellar-masslifetime relation given by the PARSEC stellar evolution model version 1.2s (Bressan et al. 2012), the theoretical stellar maximum mass limit is the most massive star that could exist at the age of the star cluster.We also plot the m ini max value predicted by the m ini max -M ecl relation ( thick dotted curve, using the GalIMF code, Yan et al. 2017).We note that the initial total stellar mass for the modeled embedded star cluster, M ecl , is transformed to M cl , the present-day mass above a completeness mass limit of m lower = 0.3 M ⊙ (the mean completeness-limited mass in our cluster sample), to compare with our observations.The large scatter in the m lower values in target clusters (Table 1) introduces a large scatter of the data points in Figure 5 (a).In addition, the transformation from M ecl to M cl accounts only for the stellar evolution, which plays the predominant role during the phase of early mass loss in a cluster (Baumgardt & Makino 2003).The combined constraints are illustrated by the thin solid curves.This indicates that, in addition to the age of a star cluster, the total cluster mass also exerts an influential role in determining a stellar mass limit following the empirical tight m ini max -M ecl relation.

On dynamical evolution
In Figure 6, we display the distribution of the powerlaw index α of the PDMF along the cluster age (left panels); and along cluster total mass (right panels).We overplot the predicted α = 2.04 (grey dashed-dotted line) from the single power-law fit of the "Kroupa IMF" (Kroupa 2001) in the mass range from 0.3 M ⊙ to 2 M ⊙ .The value of α remains stable at α ≈ 2 for clusters younger than 200 Myr, followed by a subsequent decline with increasing age (Spearman's rank correlation coefficient s = −0.16).This negative trend between α and cluster age becomes more pronounced among the population of stars within the half-mass radius r h (blue dots in panel (c)), with s = −0.40.The half-mass radius r h is defined here as the radius of the sphere that contains half of the observed total mass of the cluster above completeness limit.However, this trend is notably absent for the stars outside r h .The deviation of the value of α from the "Kroupa IMF", that is, the declining trend of α with age for stars within r h , can be attributed to the dynamical evolution of the cluster (Baumgardt & Makino 2003).As clusters become older than 100 − 200 Myr, which is a typical relaxation timescale of open clusters, an increasing number of low-mass stars escape through processes such as two-body relaxation, or being pulled away by the Galactic tidal field.Younger clusters have a higher α (steeper; a higher fraction of low-mass stars).As clusters age and approach dissolution, they experience a depletion in low-mass stars, resulting in a lower α (shallower; indicating a smaller fraction of low-mass stars).Nonetheless, the completeness-limited mass imposed by Gaia observations contributes, albeit partially, to the variability observed in the α values across clusters.Clusters with m lower below the mean value of 0.3 M ⊙ are all younger than 1 Gyr, exhibiting a mean α of 1.80 that is consistent with the initial slope of the "Kroupa IMF" between 0.2 to 2 M ⊙ with α ini 0.2−2 = 1.83.On the other hand, clusters with m lower exceeding 0.3 M ⊙ are older on average, reaching 3 Gyr in our sample.These older clusters have a mean α of 1.75 which is significantly lower than the slope of the "Kroupa IMF" in a similar mass range, being α ini 0.4−2 = 2.21 from 0.4 to 2 M ⊙ .This confirms that the PDMF of older clusters has dynamically evolved and departed more significantly from the IMF slope than that of younger clusters.
To determine the robustness of the declining trend of α with age, we derive the value of α for the PDMF of stars between the completeness limited mass and 1 M ⊙ .The correlation is identical to that of the PDMF in the left-hand panels of Figure 6.This again confirms that dynamical evolution mostly affects low-mass stars, while the dynamics of high-mass stars is governed by rapid stellar evolution.The value of α decreases when relaxation becomes dominant (see Figure 9).We do not obtain PDMF for stars above 1 M ⊙ due to low-number statistics.
A modest correlation between α and the total mass of clusters is observed.The power-law index α shows a slight increase toward more massive clusters, indicating that more massive clusters have a relatively larger fraction of low-mass stars within the selected mass range (completeness limit up to 2 M ⊙ ).This positive correlation becomes slightly more significant for stars outside r h (orange dots), a Spearman's rank correlation coefficient changing from s = 0.21 to s = 0.28.However, for stars inside r h , no such dependence is observed (s = 0.11).
The dependence of α on total cluster mass indicates that more massive clusters contain a higher fraction of faint stars.In particular, most faint stars reside in the outskirts of the massive clusters.As the stellar population evolves dynamically, low-mass stars gain kinetic energy through two-body relaxation, causing them to migrate away from the cluster center, and inhabit the outer regions.Consequently, there are more low-mass stars outside the half-mass radius r h .In low-mass clusters, the velocities of faint stars at the outskirt easily exceed the escape velocity of the cluster; these therefore quickly escape.Conversely, in more massive clusters, the deeper gravitational potential well can prevent the energetic low-mass stars in the periphery from escaping.Most young filamentary or fractal stellar groups are relatively low-mass.They are known to be in a state of disruption after gas expulsion and now are undergoing expansion (Pang et al. 2022b).The dissolution state also speeds up the evaporation of low-mass stars in these young and low-mass stellar groups.This contributes to the overall observed trend in Figure 6 (d).Pang et al. (2022b) categorized their 85 star clusters into four distinct morphological types: filamentary, fractal, halo, and tidal-tail.These morphological classifications represent a continuum of stellar density and cluster age (Pang et al. 2022b).The distribution of α values for each cluster type is summarized in Figure 7.It is observed that α values tend to decrease from filamentary to tidal-tail clusters, exhibiting a considerable amount of scatters.Halo and tidal-tail clusters are generally older (> 100 Myr) compared to filamentary and fractal type clusters (< 100 Myr), aligning with the trend seen in the dependence of α on cluster age in Figure 6 (left panels).There are two outliers in the group of filamen- In panels (c) and (d), we separate each cluster into two parts and derive the PDMF individually: within the half-mass radius (r < r h : blue dots) and outside half-mass radius (r > r h : orange dots).The orange and blue triangles are computed in the same manner as in panels (a) and (b).We exclude the disrupted cluster Group X from panels (c) and (d) since it is hard to define the cluster center for its two-piece-fragmented shape.The quantity s is Spearman's rank correlation coefficient, and p is the probability of the null hypothesis (i.e., that no correlation exists between two variables) of the correlation test.A p value of less than 0.1 indicates that the null hypothesis is rejected.m lower and m upper = 2 M ⊙ , and therefore have larger uncertainties in α.

On metallicity
In this section, we explore the correlation between the PDMF and the metallicity ([M/H]) of star clusters.The metallicities for our sample stars are taken from Li et al. (2023b), who utilized a neural-network-based regression model AspGap to process Gaia DR 3 XP spectra.This model was trained using APOGEE spectra data (Abdurro'uf et al. 2022).We also cross-match the member stars in our target clusters with the data release 17 of APOGEE survey (Abdurro'uf et al. 2022) to obtain their chemical abundances.To ensure robust statistics, we consider only the clusters with a minimum of three stars that have measurable metallicity, allowing us to compute both the mean and dispersion of metallicity within each cluster.
To estimate the mean and standard deviation of the metallicity ([M/H]) for each cluster, we apply the Maximum Likelihood Estimation (MLE) method.We assume that the distribution of [M/H] values follows a Gaussian distribution for each cluster.The MLE technique enables us to identify the most probable mean (µ [M/H] ) and standard deviation (σ [M/H] ) values for this distribution based on the observed data.
In constructing the likelihood function, we take into account the number of member stars in each cluster (N ), the measured [M/H] values (x i ), and the individual measurement uncertainty (δ i ).By maximizing this function, we can determine the values of the parameters that provide the optimal fit to the observed metallicity.We compute the mean [M/H] value and the inherent dispersion that maximizes this likelihood function: .
(5) We note that the metallicity measurements are subject to uncertainties, particularly for low-temperature stars.Commonly used stellar atmosphere models are optimized for FGK stars.Obtaining atmospheric properties of M-type stars is complex due to the presence of numerous broad absorption features from molecular lines that alter the atmosphere's structure (Iyer et al. 2023).This issue extends to the MARCS models utilized in the APOGEE survey.A distinct pattern is evident in the [M/H]-T eff diagram for the stars with T eff < 4500 K. Therefore, we exclude stars with T eff < 4500 K when establishing the relationship between the power-law index α and metallicity, to avoid large uncertainties.The measured values of [M/H] and the corresponding dispersions for the target clusters are presented in columns 7 and 8 in Table 1.
In Figure 8, the dependencies of the PDMFs on metallicity obtained from stars with effective temperature T eff > 4500 K are presented.No discernible trend between metallicity and α is seen in both Gaia DR 3 data and APOGEE data.In contrast to both theoretical and observational investigations (Skillman 2008;Martín-Navarro et al. 2015;Smith 2020;Yan et al. 2020;Sharda & Krumholz 2022;Li et al. 2023) that propose a strong correlation between metallicity and the low-mass IMF slope, our sample does not provide evidence for such a relationship.This discrepancy may arise from the predominant similarity in metallicity across most of the clusters in our sample, resulting in negligible variations in the IMF.In particular, the PDMF may not be sensitive to the expected low-mass IMF differences that might be washed out by the pronounced effects of dynamical evaporation.
An additional noteworthy pattern emerges, indicating that younger clusters (blue to magenta dots) exhibit lower metallicity compared to their older counterparts (orange to yellow dots).The older and more enriched clusters potentially might not have originated in the solar neighborhood; rather, they might have migrated from the inner disk to their current local posi-tions through radial migration processes (Minchev et al. 2009).
Recently, Pang et al. (2022c) discovered the Collinder 132-Gulliver 21 Stream (270 pc long) in the solar neighborhood, with stellar populations showing an age difference up to 250 Myr.The oldest generation in the stream is approximately 0.3 dex (based on Gaia DR3 data) more metal-rich than the youngest generation.This may provide additional evidence of radial migration in the solar neighborhood.A similar pattern is also found in other open clusters at larger distances from the Sun (see, e.g., Magrini et al. 2023).4.3.On black holes Torniamenti et al. (2023) proposed that stellar-mass black holes (BHs) can be retained even in open clusters with low escape velocities (< 3 km s −1 ).Their models incorporated 2-3 BHs in the Hyades cluster, leading to an increased central velocity dispersion of observed star clusters.However, distinguishing this BH signature from the effects of binary stars remains challenging.
We do not observe stars that are sufficiently massive to be candidate BH progenitors.Given that most clusters in our sample are older than 10 Myr, BH progenitors should already have evolved.The most massive stars in the young (< 10 Myr) star clusters, on the other hand, are not massive enough to evolve into BHs.However, this cannot exclude the that BHs may have been present in the sample of young clusters in the past, because the massive BH progenitors can be ejected due to dynamical 3-body and 4-body interactions (Oh et al. 2015).Even for the clusters that lack BH progenitors from stellar evolution, there is still the possibility to produce stellarmass BH via binary mergers and frequent physical collisions (Oh & Kroupa 2018;Spera et al. 2019;Di Carlo et al. 2020).Consequently, we suggest at most a few BHs may exist in our sample of open clusters.
In globular clusters, on the other hand, stellar-mass BHs form a subsystems (Breen & Heggie 2013).When a BH subsystem contains more than 40 BHs, this system can significantly impact the cluster dynamics.Baumgardt et al. (2023) carried out star cluster simulations with retention of a BH subsystem by adopting lowvelocity kicks.Their findings reveal that the presence of a BH subsystem significantly accelerates the two-body relaxation between BHs and massive stars.During the phase of fast relaxation, approximately 15 to 20 percent of the most massive stars are lost.Consequently, the retention of BHs in star clusters primarily amplifies the escape rate of higher-mass stars, and this effect can be observed through PDMF.
We compare the power-law index α of the PDMF from our cluster sample with the models of Baumgardt et al. (2023) in Figure 9.The models (plus signs) align with the lower edge of our cluster sample.Differentiation between various BH retention models primarily occurs in dynamically old clusters (age > relaxation time).The relaxation time of each cluster is computed with equation (7.108) in Binney & Tremaine (2008), using the half-mass radius, the total mass, and the number of members in each cluster.However, the majority of clusters in our sample are younger than one relaxation time, indicating a scarcity of old open clusters in the solar neighborhood for model comparison.
The solar neighborhood, spanning a radius of 500 pc, resides within the Local Arm (Reid et al. 2019) and is characterized by abundant young star formation regions and young star clusters.The older clusters in our sample likely did not form in situ, but instead underwent radial migration from the inner disk (as discussed in Section 4.2).Given that the typical survival timescale of open clusters is on the order of a hundred million years, the observed aged open clusters are skewed toward those that were initially massive at birth or were located in an orbit that experienced less frequent tidal interactions, such as those positioned in the direction of the Galactic anti-center.That is, old clusters are mostly to be found at larger distances.However, the PDMF of distant old open clusters is significantly influenced by the completeness mass limit of Gaia data.This should be properly addressed to facilitate meaningful comparisons with nearby clusters in the sample.
Note that in Baumgardt et al. (2023), a BH population comprises a few dozen BHs.The presence of this many BHs can significantly influence the evolution of the system.On the other hand, in open clusters (which typically only contain a few black holes, Torniamenti et al. 2023), the influence of these stellar-mass BHs on the PDMF is limited.Instead, the impact of these BHs is similar to that of massive compact binary systems in the core: relaxation due to gravitational interactions with neighbouring stars.Additionally, due to the higher metallicity of open clusters compared to globular clusters, the masses of stellar-mass BHs formed in open clusters are considerably lower (Belczynski et al. 2010).This further weakens the influence of BHs in open clusters.Finally, the influence of a stellar-mass BH on the stellar kinematics only becomes apparent on a relatively long time scale.The older clusters are therefore preferred targets to investigate the potential presence of BH candidates.Thus, distant, old, and massive open clusters could provide an opportunity to detect stellar-mass BHs.We analyze the PDMF of 93 star clusters using Gaia DR 3 data.The stellar membership of 92 of the clusters is taken from previous studies (Pang et al. 2021a(Pang et al. ,b, 2022b,c;,c;Li et al. 2021).These studies determined stellar membership using the machine learning algorithm StarGo.Membership of the remaining cluster, NGC 6991, is identified in this work.We evaluate the impact of unresolved binary systems on stellar mass estimation.We adopt three mass ratio profiles: a uniform profile, the K07 profile (Kouwenhoven et al. 2007), and the F05 profile (Fisher et al. 2005) for the correction of the observed stellar mass.A binarity probability is assigned to observed stars, based on the empirical relation between the binary frequency and primary mass taken from Offner et al. (2022).
We characterize the PDMF with a power-law index, α, which is determined using a maximum likelihood method that avoids biases introduced by data binning.The likelihood function is formulated for each cluster, allowing for robust statistical analysis.The value of α is only fitted for a given mass range: between the Gaia completeness limited mass (lower limit) and a consistent upper mass limit of 2 M ⊙ for all clusters.The resulting histograms of α values, considering different mass ratio distributions for binary corrections, exhibit negligible differences.
We explore the relationship between the maximum observed stellar mass, m max , and the observed total mass above the completeness limit of the star clusters.The analysis reveals a positive correlation between m max and the total cluster mass, where more massive young clusters tend to host stars with a higher m max .Older clusters currently have smaller m max values due to stellar evolution.Our results suggest that, beyond the age of star clusters, their total mass plays a significant role in determining the maximum stellar mass limit, reinforcing previous findings in the literature (Yan et al. 2023b).
Based on the PDMF of these 93 star clusters, we investigate the dependence of the power-law index, α, on dynamical evolution, metallicity, and black holes.
• There is a notable correlation between the powerlaw index, α, of the PDMF and the age of the clusters.Specifically, the α value for clusters younger than 200 Myr remains relatively stable at a value of α ≈ 2, which is consistent with the α value for a single power-law fit of the canonical Kroupa (2001) IMF in the mass range between 0.3 and 2 M ⊙ , but declines for older clusters.This negative trend becomes more pronounced when considering stars within the half-mass radius.The dynamical processes within the cluster, such as two-body relaxation and the Galactic tidal stripping, contribute to the removal of low-mass stars and produce the observed variation in α.
• The dependence of the PDMF slope on the total cluster mass is weak, with a slight increase in the α value for more massive clusters.This dependence becomes slightly more pronounced when focusing on stars outside the half-mass radius.As low-mass stars migrate to the outer region of the cluster via two-body relaxation, the cluster potential determines whether they will remain bound to the cluster or whether they will escape.Less massive clusters exhibit a higher rate of escaping faint stars as a result of a lower escape velocity.
• Regarding cluster morphology, we find a decreasing trend in the values of α from filamentary to fractal, to halo and to tidal-tail clusters, which are also a sequence in terms of increasing cluster age.This resembles the correlation between the α value and the cluster age and again emphasizes the importance of secular dynamical evolution in shaping the PDMF of star clusters.
• Utilizing data from Li et al. (2023b) and APOGEE survey, we examine the relation between the metallicity ([M/H]) of cluster members and the power-law index α of the PDMF.We exclude stars with T eff < 4500K K due to uncertainties in metallicity measurements for these stars.Our results reveal no discernible trend between the value of α and the mean metallicity derived from stars with T eff > 4500K.However, a distinct pattern emerges: younger clusters exhibit a lower metallicity compared to their older counterparts.This discrepancy is attributed to potential radial migration, indicating that the older, more metal-rich clusters may have originated from the inner circle of the solar neighborhood and subsequently migrated outward.
• We attempt to explore the role of BHs in shaping star cluster PDMF, building on the findings of recent studies.We study the evolution of the α value with age (measured in relaxation times), and make a comparison between our observational results with the numerical simulations by Baumgardt et al. ( 2023) that included a BH population in their star cluster models.The differentiation between BH retention models is more pronounced in dynamically old clusters, whereas the majority of the analyzed star clusters are younger than one relaxation time.Therefore, the current star cluster sample cannot provide strong constraints for these models.There is a need for more distant, older, and more massive open clusters to explore the potential presence and impact of stellar mass BHs.
Our study provides valuable insights into the PDMF of star clusters in the solar neighborhood.The PDMF is heavily shaped by the internal and external dynamical evolution.The predominantly younger age range of our cluster sample emphasizes the need for observations of older and more massive open clusters that are located further away from the solar neighborhood, to refine star cluster models that include stellar mass black holes.The incompleteness limit from observations will need to be carefully considered when investigating distant old clusters.

Figure 1 .
Figure1.Distributions of mass ratio q for binary systems.All three mass ratio distributions are normalized to unity.The uniform distribution (Uni) is shown as the blue line,Kouwenhoven et al. (2007, K07)  is shown as the red dashed curve, andFisher et al. (2005, F05)  as the orange dotted curve.

Figure 2 .
Figure 2. Dependence of binary probability p on primary mass m ′ .The red crosses are empirical values obtained from Offner et al. (2022).The orange curve is the function that is interpolated by the red crosses.

Figure 3 .
Figure3.The most probable PDMFs of example clusters Pleiades (a, b, c) and Praesepe (d, e, f) after binary correction, considering three different mass ratio distributions.The red vertical dashed line indicates the completeness limited mass of Gaia DR 3 data m lower , which is 0.28 M⊙ and 0.31 M⊙ for Pleiades and Praesepe respectively.The red vertical dotted line corresponds to the upper mass limit mupper of 2 M⊙.The blue curves are the mass distributions of Pleiades and Praesepe.The PDMFs (orange dashed lines) are determined only for the stellar mass between m lower and mupper.The power-law index α value obtained from the maximum likelihood method in Section 3.2 is indicated in each panel.

Figure 4 .
Figure 4. Histograms of the power-law indices (α) of the PDMFs for a total of 93 star clusters.Histograms are shown for results obtained after binary corrections with different mass ratio distributions.The blue histogram shows the results obtained for a uniform (Uni) mass ratio distribution; the red dashed histogram for the mass ratio profile ofKouwenhoven et al. (2007, K07); and the orange dotted histogram for that ofFisher et al. (2005, F05).

Figure 5 .
Figure 5. (a): Relation between mmax, the maximum stellar mass in each cluster, and M cl , the total observed cluster mass above the completeness limit.Each dot represents one open cluster, and its color represents the cluster age, which is indicated by the color bar.The grey triangles in panel are average values of the total cluster mass and the maximum stellar mass for all 15 or 16 clusters in each bin, with the standard deviation indicated by the error bar.The horizontal grey dashed-dotted line indicates the upper mass limit, mupper = 2 M⊙, adopted in the PDMF derivation.The thin dashed horizontal lines present the theoretical stellar maximum mass limit at a given age.The thick dotted curves indicate the values of m ini max for embedded star clusters given by the relation m ini max -M ecl (Yan et al. 2017).The combined constraints are illustrated by the thin solid curves.(b): Relation between mmax and the cluster age.The total cluster mass above the completeness limit, M cl , is indicated by the color bar.

rFigure 6 .
Figure6.Dependence of the power-law index (α) of the most probable PDMF on the cluster age (left panels) and total cluster mass above completeness limit, M cl (right panels).The uniform mass ratio distribution is used for binary correction.The grey dashed-dotted line in each panel corresponds to α = 2.04 computed from theKroupa (2001)  IMF in the mass range from 0.3 M⊙ to 2 M⊙.The values of α in panels (a) and (b) are computed for all members in the cluster.The orange triangles are average values of cluster age and α for all 15-16 clusters in each bin, with the standard deviation indicated by the error bar.In panels (c) and (d), we separate each cluster into two parts and derive the PDMF individually: within the half-mass radius (r < r h : blue dots) and outside half-mass radius (r > r h : orange dots).The orange and blue triangles are computed in the same manner as in panels (a) and (b).We exclude the disrupted cluster Group X from panels (c) and (d) since it is hard to define the cluster center for its two-piece-fragmented shape.The quantity s is Spearman's rank correlation coefficient, and p is the probability of the null hypothesis (i.e., that no correlation exists between two variables) of the correlation test.A p value of less than 0.1 indicates that the null hypothesis is rejected.

Figure 7 .
Figure 7. Dependence of the power-law index (α) of the most probable PDMF on cluster morphologies: f1: filamentary (blue), f2: fractal (red), h: halo (purple), and t: tidal tail (green).The colored rectangles indicate the inner quartile range (IQR, 75th percentile minus 25th percentile).The median value is indicated with a horizontal line.The upper and lower whiskers of each colored rectangle mark the maximum and minimum values within 1.5 IQR.Clusters outside 1.5 IQR (whiskers) are outliers indicated as open circles: Stock 23 and UBC 31 gp1 (blue circles).

Figure 8 .Figure 9 .
Figure 8. Distributions of the power-law index (α) of the most probable PDMF and cluster metallicity [M/H].The age of the cluster is indicated by the color bar on the right.The metallicity of clusters [M/H] in (a) and (b) is obtained using members that have [M/H] available from Gaia DR 3 (a) and APOGEE (b), and that also have T eff > 4500 K. Faint stars with T eff < 4500 K have been excluded due to their large uncertainties.

Table 1 .
PDMF and related parameters for 93 star clusters.

Table 1 continued
Table 1 (continued) fitting.The age of NGC 6991 is obtained in this work (see Appendix A