IMF-aware N-body forward modelling of membership selection functions for young open clusters

Introduction

The stellar initial mass function (IMF), the distribution of stellar masses at birth, sets the luminosity, lifetime and feedback of a stellar population, and is one of the largest sources of systematic uncertainty in stellar evolution models (Salpeter 1955; Hoppins 2018; Kroupa et al. 2026). It has been described by different mathematical forms, such as multi-part or tapered power laws (Kroupa 2001; De Marchi et al. 2010) and log-normal distributions (Miller & Scalo 1979; Chabrier 2003). Beyond the choice of mathematical form, a long-standing question is whether the IMF is universal (Bastian et al. 2010) or varies systematically with star-forming density and metallicity (Marks et al. 2012; Li et al. 2023). Star clusters are well suited to testing these possibilities, because their members share a common age, distance and chemical composition, and clusters from different environments can be compared (Lada & Lada 2003; Portegies Zwart et al. 2010; Krumholz et al. 2019). Young nearby open clusters are particularly useful: when their ages are short compared with their relaxation times, mass segregation and the preferential loss of low-mass stars have had less time to alter their mass functions (Spitzer 1987; Baumgardt & Makino 2003), and Gaia astrometry [20] allows their mass functions to be measured over a broad mass range (e.g. Cordoni et al. 2023). Recent $N$-body modelling that accounts for dynamical evolution and unresolved bi-

© 2026 The Authors naries suggests tentative cluster-to-cluster scatter in the high-mass IMF slope (Hobart et al. 2026). Whether this scatter reflects genuine IMF variations also depends on how well the incompleteness and field star contamination introduced by membership selection are quantified.

We therefore focus on this selection effect, which affects the mass function derived from any star counting approach. Incompleteness and contamination that are independent of magnitude would only rescale the mass function, whereas those that vary with magnitude, and hence with stellar mass, distort its shape. Faint stars have larger astrometric uncertainties (Lindegren et al. 2021), which broaden their observed distributions in proper motion and parallax. In addition, low mass stars are preferentially found in the outskirts of mass segregated clusters (Tarricq et al. 2022), and many young clusters are observed to be expanding (Della Croce et al. 2024), reducing the density contrast of their members against the field. Both effects make faint members harder to separate from field stars. Correcting the mass function of a cluster therefore requires the completeness $f_{\mathrm{comp}}(G)$, the fraction of true members that are recovered, and the contamination $f_{\mathrm{cont}}(G)$, the fraction of selected stars that are field stars, as functions of magnitude. We refer to these as the membership selection functions; they are specific to both the cluster and the procedure used to identify its members. They differ from the selection function of the Gaia catalogue itself, which describes whether a star is included in the catalogue (Cantat-Gaudin et al. 2023).

Gaia has enabled homogeneous member lists for thousands of open clusters (Cantat-Gaudin et al. 2018; Castro-Ginard et al. 2020, 2022; Hunt & Reffert 2023; Perren et al. 2023), constructed with a range of membership algorithms [8] (Hunt & Reffert 2021; Pera et al. 2021; Van Groeningen et al. 2023). The membership probabilities provided with these lists are specific to each algorithm and are not, in general, calibrated estimates of the completeness and contamination of the lists. Hunt & Reffert (2024) estimated magnitude-dependent completeness for individual clusters and used it to correct their masses, but modelled member recovery from astrometric uncertainties alone, without field star contamination. Hunt et al. (2026) went further: in addition to the probability that a cluster as a whole is detected, they illustrated both completeness and contamination at the membership level with injection–recovery tests, but their synthetic clusters had a fixed Kroupa IMF (Kroupa 2001) and King spatial profiles (King 1962), without dynamical evolution.

Dedicated mass function studies have typically mitigated these effects rather than directly correcting the selection functions of individual member stars. For example, some studies restrict their analysis to magnitude ranges where the survey is complete (Olivares et al. 2018). Others select members conservatively and verify their completeness against deeper surveys (Hobart et al. 2026), or use simulated clusters to demonstrate that their member lists are sufficiently complete and pure without correction (Jiang et al. 2024). Alternatively, some studies avoid strict membership selection altogether by fitting field stars as a separate component in the colour–magnitude diagram (Li et al. 2025, 2026). For existing member stars, completeness and contamination as functions of magnitude, including the effects of cluster evolution and structure, therefore remain largely unquantified.

Measuring these selection functions requires the true membership of every star, which is unknown for observed clusters but known by construction in simulations. Simulating a cluster, however, requires assumptions about its stellar population and structure. The IMF is particularly problematic: if the selection functions depend on it, corrections based on a fixed IMF would partly assume the very quantity they are meant to constrain. The IMF can affect the selection functions directly, by setting how many members lie at each magnitude relative to field stars, and indirectly through the evolution of the cluster, since an IMF richer in massive stars leads to stronger stellar evolution mass loss and expansion [36] (Chernoff & Weinberg 1990; Haghi et al. 2020; Wang & Jerabkova 2021). Capturing this indirect effect requires evolving each simulated cluster from an initial state, rather than imposing a fixed present-day profile. In addition, the selection functions depend on the cluster, its surrounding field and the membership procedure (Hunt & Reffert 2021, 2024), so those derived from generic simulated clusters cannot be transferred directly to a particular member list.

$N$-body models of individual clusters have been used to infer their initial conditions and IMF (Baumgardt et al. 2023; Weis et al. 2025; Hobart et al. 2026) and to calibrate their initial masses (Zheng et al. 2026), but the simulated stars are compared with the observations without passing through the same membership selection as the observed members. Deriving selection functions for a given member stars therefore requires cluster-specific $N$-body models with the IMF as an explicit input, whose simulated stars pass through the same membership selection as the observed members.

We construct such models with an IMF-aware $N$-body forward-modelling framework. For each cluster and each assumed IMF, we evolve an $N$-body model to the cluster age, observe it with the Gaia uncertainties of that cluster and a mock field population, and identify its members with the same clustering algorithm as used for the observed list. The initial cluster mass is calibrated so that the number of recovered members matches the observed one, which fixes the initial mass for each assumed IMF. The completeness and contamination are then measured from the known origin of each simulated star, with their scatter estimated from independent realizations. The framework can be applied to member lists from other algorithms by replacing only the membership identification step. We apply it to 35 nearby young open clusters with HDBSCAN member lists (Hunt & Reffert 2023, 2024), using three IMF slopes above $0.5\,M_{\odot}$ ($\alpha_2=-2.7$, $-2.3$, and $-1.9$), and examine how the selection functions vary with magnitude and between clusters, whether they depend on the assumed IMF, whether corrections derived with a mismatched IMF bias the recovered luminosity function, and how these differences relate to the evolution and present-day structure of the models.

The paper is organized as follows. § 2 describes the observational sample to which we apply the framework. § 3 introduces the IMF-aware $N$-body forward-modelling framework and details its construction. § 4 presents the magnitude-dependent and global completeness and contamination for different IMF slopes, their effect on the recovered present-day luminosity function, and their dependence on cluster properties. § 5 discusses the evolutionary and structural changes associated with the IMF dependence, the dependence on the membership algorithm, and the caveats of this work. § 6 summarizes our conclusions.

Observational sample

Our observational sample consists of 35 nearby young open clusters with distances $\lesssim600\mathrm{pc}$ and ages $\lesssim100\mathrm{Myr}$. The targets are initially selected by cross-matching the clusters studied by Pang et al. (2022) with the Gaia DR3 catalogue of Hunt & Reffert (2024). Of the 37 matched objects, Alessi24 and BH164 are excluded because the mock-field data required for their simulations were unavailable. All observational member lists and reference cluster properties adopted in this work are taken from Hunt & Reffert (2023, 2024). The identifications and adopted properties of the 35 clusters are listed in Table A1 and described in Appendix A.

HR24 name	P22 name	R.A. (deg)	Decl. (deg)	$d$ (pc)	$\mu_{\alpha^\ast}$ (mas yr$^{-1}$)	$\mu_{\delta}$ (mas yr$^{-1}$)	$t$ (Myr)	$M_{\mathrm{cl}}$ ($M_\odot$)	$N_{\mathrm{cat}}^{\mathrm{obs}}$
ASCC 105	ASCC 105	295.472	27.386	553.4	1.370	-1.618	$93.4^{+73.0}_{-40.9}$	$307.3 \pm22.4$	173
ASCC 127	ASCC 127	347.302	65.183	374.2	7.410	-1.645	$12.2^{+7.1}_{-5.1}$	$168.6 \pm20.4$	247
Briceno 1	ASCC 16	81.083	1.696	342.1	1.384	-0.016	$8.2^{+3.1}_{-3.8}$	$118.5 \pm19.3$	171
ASCC 19	ASCC 19	81.978	-1.856	352.1	1.168	-1.182	$7.4^{+3.6}_{-2.7}$	$56.1 \pm13.0$	61
ASCC 58	ASCC 58	153.706	-55.020	471.4	-13.302	2.760	$57.0^{+48.1}_{-26.7}$	$266.5 \pm32.9$	236
Alessi 20	Alessi 20	2.643	58.757	424.3	8.117	-2.427	$6.7^{+2.3}_{-2.8}$	$204.7 \pm26.6$	252
Alessi 5	Alessi 5	160.902	-61.112	391.8	-15.400	2.544	$56.4^{+40.1}_{-22.7}$	$657.9 \pm53.6$	571
BH 99	BH 99	159.554	-59.109	440.2	-14.477	0.980	$44.8^{+42.7}_{-19.1}$	$739.4 \pm60.2$	643
Collinder 135	Collinder 135	109.406	-36.921	295.9	-10.078	6.230	$29.8^{+23.9}_{-13.1}$	$128.3 \pm25.7$	209
Collinder 140	Collinder 140	110.802	-32.012	376.2	-8.088	4.756	$17.3^{+12.6}_{-7.7}$	$164.2 \pm22.8$	203
Collinder 69	Collinder 69	83.814	9.888	393.8	1.215	-2.047	$5.2^{+1.8}_{-1.8}$	$510.6 \pm92.7$	1247
Theia 13	Gulliver 6	83.132	-1.753	408.8	0.104	-0.222	$6.4^{+3.1}_{-2.4}$	$202.3 \pm25.9$	249
IC 2391	IC 2391	130.260	-53.042	150.2	-24.797	23.267	$27.3^{+20.2}_{-12.9}$	$203.9 \pm27.6$	376
IC 2602	IC 2602	160.973	-64.393	150.6	-17.775	10.669	$26.0^{+21.6}_{-10.8}$	$344.2 \pm42.3$	638
IC 4665	IC 4665	266.532	5.635	343.3	-0.881	-8.532	$32.7^{+20.2}_{-14.4}$	$278.9 \pm25.2$	308
LP 2383	LP 2383	95.380	-16.097	365.0	-5.370	5.049	$37.9^{+40.9}_{-19.6}$	$322.4 \pm33.9$	422
ZHBJZ 1	LP 2439	103.947	-5.791	283.7	-7.318	-2.445	$13.5^{+8.8}_{-6.9}$	$275.0 \pm37.2$	490
UPK 640	LP 2442	249.692	-39.507	173.9	-11.882	-21.391	$9.6^{+4.8}_{-4.4}$	$484.1 \pm44.0$	764
NGC 1980	NGC 1980	83.801	-5.921	378.3	1.265	0.472	$6.9^{+2.9}_{-2.6}$	$245.9 \pm35.5$	364
NGC 2232	NGC 2232	96.978	-4.772	318.0	-4.717	-1.770	$13.3^{+9.2}_{-5.6}$	$185.8 \pm25.5$	286
NGC 2451A	NGC 2451A	115.971	-38.241	190.1	-20.987	15.327	$26.1^{+22.7}_{-12.9}$	$220.8 \pm30.4$	407
NGC 2451B	NGC 2451B	116.176	-37.942	367.5	-9.400	4.744	$27.6^{+19.5}_{-10.7}$	$494.1 \pm52.0$	615
NGC 2547	NGC 2547	122.533	-49.195	382.4	-8.568	4.331	$21.7^{+19.5}_{-9.2}$	$530.4 \pm54.3$	675
NGC 3228	NGC 3228	155.278	-51.845	477.1	-14.856	-0.595	$30.3^{+25.8}_{-12.1}$	$250.8 \pm23.2$	200
NGC 6405	NGC 6405	265.067	-32.228	451.9	-1.368	-5.847	$52.4^{+38.5}_{-22.3}$	$1531.7 \pm173.4$	974
NGC 7058	NGC 7058	320.409	50.830	362.8	7.429	2.805	$69.7^{+52.9}_{-35.3}$	$334.1 \pm43.2$	424
NGC 7429$^{\mathrm{m}}$	RSG 7	344.090	59.912	417.4	4.771	-2.245	$78.1^{+58.5}_{-32.0}$	$71.8 \pm14.8$	67
RSG 8	RSG 8	345.453	59.509	481.4	5.664	-1.716	$13.6^{+8.8}_{-5.5}$	$175.1 \pm19.1$	151
Roslund 5	Roslund 5	302.617	33.741	538.9	2.052	-1.236	$110.1^{+95.0}_{-41.8}$	$489.1 \pm43.1$	360
Stephenson 1	Stephenson 1	283.578	36.881	354.7	1.121	-3.018	$26.6^{+23.5}_{-12.1}$	$267.7 \pm27.2$	322
Stock 23	Stock 23	49.184	60.071	598.2	-4.324	-0.993	$71.8^{+46.7}_{-30.9}$	$203.3 \pm37.4$	106
Alessi-Teutsch 10	UBC 19	56.331	29.742	392.6	2.671	-5.238	$3.4^{+0.5}_{-0.9}$	$46.5 \pm9.5$	85
OC 0280$^{\mathrm{m}}$	UBC 31	60.521	32.178	368.0	3.628	-5.259	$7.8^{+3.6}_{-3.3}$	$52.4 \pm12.5$	78
Alessi 36	UBC 7	106.565	-37.596	275.3	-9.669	7.046	$28.9^{+31.9}_{-14.3}$	$105.8 \pm16.8$	198
UPK 82	UPK 82	298.464	26.202	537.9	2.278	-2.169	$77.9^{+72.3}_{-35.1}$	$194.7 \pm26.5$	142

This work focuses on the demonstration of the method that can provide the IMF-aware model for the membership selection function, which plays an important role in mass function determination of a given star cluster, rather than on handling as many star clusters as possible. Therefore, the sample needs to contain sufficient clusters to demonstrate the method, and should also consist predominantly of young clusters, which are critical for future IMF measurement studies.

The observational member lists adopted in this work were obtained with HDBSCAN in the all-sky search of Hunt & Reffert (2023). HDBSCAN² (Hierarchical DBSCAN, Campello et al. [8]) combines the density-based approach of DBSCAN (Density-Based Spatial Clustering of Applications with Noise, Ester et al. (1996)) with hierarchical clustering to effectively handle datasets of varying densities. To search for clusters closer than 750 pc, Hunt & Reffert (2023) divided the sky into 192 level-2 HEALPix (Hierarchical Equal Area isoLatitude Pixelation) pixels and selected sources with $\varpi> 1$ mas. Each pixel and its eight nearest neighbours formed a field of approximately $44^{\circ} \times44^{\circ}$, with adjacent fields overlapping. For clusters closer than 150 pc, they instead used a single all-sky field containing stars with photo-geometric distances below 250 pc (see their Section 2.3). Within each HEALPix field, positions, proper motions, and parallaxes were transformed to a coordinate frame centred on the field and rescaled so that all five dimensions had the same median and interquartile range. For the all-sky field, positions and parallaxes were converted to Cartesian coordinates using photo-geometric distances, and proper motions to tangential velocities (their Section 3.1). They adopted the same HDBSCAN hyperparameter values as in their previous study (Hunt & Reffert (2021)), with $m_{\mathrm{clSize}} \in\{10, 20, 40, 80\}$ and $m_{\mathrm{Pts}} = 10$. Detections were required to exceed $3\sigma$ in the cluster significance test before the results from different hyperparameter runs and fields were combined and duplicate detections removed. During this merging, only clusters closer than 175 pc were retained from the all-sky search, whereas those from the level-2 HEALPix search were retained only between 150 and 750 pc. The overlap allowed the better detection to be selected for clusters near the boundary (see their Section 3.2). The minimum cluster size $m_{\mathrm{clSize}}$ adopted for each object is recorded in the online cluster catalogue, as described in their Table A.1. We use the same HDBSCAN algorithm to recover members in our forward models, as described in § 3.5.1.

Method

Forward-Modelling Framework

We construct an IMF-aware forward model that follows each observed cluster from an assumed initial state to a present-day membership catalogue, as shown in Figure 1. The model combines $N$-body simulation (step 1), synthetic Gaia observations with a mock field population (step 2), and HDBSCAN membership recovery (step 3). The IMF is an explicit input: for each IMF slope $\alpha_2$, we generate a new set of stellar masses and follow the stellar and dynamical evolution of the resulting cluster in a separate $N$-body model. The output is a simulated counterpart of the observed member list, subject to the same magnitude and footprint cuts as the observed sample. As in Zheng et al. (2026), we calibrate $M_{\mathrm{cl},0}$ by matching simulated and observed star counts, but after membership recovery rather than within the tidal radius, so that both counts are subject to incompleteness and contamination. We repeat steps 1–3 with trial values of $M_{\mathrm{cl},0}$, updating $r_{\mathrm{h},0}$ through the adopted mass–radius relation (§ 3.2), until the selected simulated and observed sample counts agree to within 20 per cent (§ 3.5). In step 4, we use the known origin of each simulated star in the accepted model to measure the membership selection functions, namely completeness and contamination functions of $G$ magnitude, and repeat the calculation across independent realizations to estimate their scatter.

Vector (PDF)

For a given $\alpha_2$, the same forward model, which follows the cluster from its initial state to the present day, thus serves both to calibrate the initial mass and to derive the selection functions. For each cluster, we obtain an accepted model for each $\alpha_2$, all matching the same observed member count within the adopted tolerance. Comparing their selection functions allows us to examine whether and how they vary with the assumed IMF for that cluster.

Initial Conditions and N-body Evolution

We generate the initial cluster models with a modified version of MCLUSTER³ (Wang et al. (2019); Küpper et al. [41]). We adopt a continuous, two-part power-law IMF (Kroupa (2001)),

where $k$ is the normalization constant. We consider $\alpha_2 = -1.9$, $-2.3$, and $-2.7$ to compare membership selection functions under different IMFs. We assign stellar masses by optimal sampling (Kroupa et al. (2013)) to reduce stochastic variations in the evolution of low-mass clusters (Wang & Jerabkova (2021)). For each $M_{\mathrm{cl},0}$ and $\alpha_2$, we recalculate the discrete mass spectrum and maximum initial stellar mass $m_{\max}$. In our implementation, $k$ and $m_{\max}$ jointly satisfy

Here $150\,M_{\odot}$ is the IMF upper limit, while $m_{\max}$ is the mass of the most massive star assigned to the model.

We calibrate the initial cluster mass $M_{\mathrm{cl},0}$ separately for each $\alpha_2$ (§ 3.5) and set the initial half-mass radius using the mass–radius relation of Marks & Kroupa (2012), $r_{\mathrm{h},0} = 0.1 \times(M_{\mathrm{cl},0}/M_{\odot})^{0.13}$ pc. For our calibrated masses of $M_{\mathrm{cl},0} \approx25$–$1800\,M_{\odot}$, this gives $r_{\mathrm{h},0 \approx0.15}$–$0.26$ pc. Models adopting this relation may remain more compact at the present day than the observed clusters (Zheng et al. (2026)). We nevertheless retain this relation to limit computational cost and avoid parameter degeneracies that would be difficult to resolve if $r_{\mathrm{h},0}$ were varied independently alongside $M_{\mathrm{cl},0}$ and $\alpha_2$. The clusters are initialized as King models with $W_0 = 5$ (King (1966)), scaled to $r_{\mathrm{h},0}$ and in virial equilibrium, without primordial mass segregation, fractal substructure, or primordial binaries.

We evolve the clusters with the $N$-body code PeTar⁴ (Wang et al. (2020b)), which combines the particle–particle/particle–tree ($\mathrm{P^3T}$) method (Oshino et al. (2011)) with slow-down algorithmic regularization (SDAR; Wang et al. (2020a)) and is parallelized with the Framework for Developing Particle Simulators (FDPS; Iwasawa et al. (2016, 2020)). Stellar and binary evolution, including wind mass loss, stellar-type transitions, mass transfer, and mergers, is followed with the updated SSE/BSE codes (Hurley et al. (2000, 2002); Banerjee et al. 2020), and the Galactic potential is included through galpy (Bovy 2015).

We adopt $Z = 0.02$ for stellar evolution and MWPotential2014 for the Galactic potential. For each cluster, a centre-of-mass particle is integrated backwards from the present-day Galactocentric position and velocity derived from the catalogue of Hunt & Reffert (2024). The backward integration time is the catalogue median age rounded down to the nearest Myr (Appendix A). The $N$-body models for all mass and IMF trials of the cluster start from the resulting position and velocity and are evolved forward over this interval to the present day. In the final snapshot, we retain all particles, including remnants and escapers still in the simulation, without any radial or bound-membership cut. This snapshot is the input to the synthetic Gaia observations described in § 3.3.

Synthetic Gaia Observations and Footprint

We use GALEVNB (Pang et al. 2016) to compute absolute $G$, $G_{\mathrm{BP}}$, and $G_{\mathrm{RP}}$ magnitudes only for particles with star. type $< 10$ (i.e. main-sequence, Hertzsprung-gap, first-giant-branch, core-helium-burning, and asymptotic-giant-branch stars, and naked helium stars). We assign all stars in each cluster the catalogue median $A_V$ (AV50) from Hunt & Reffert (2024) and convert it into extinction in each band using the corrected polynomial coefficients of Monetiro et al. (2024), with extrapolation where necessary. We then obtain noiseless apparent magnitudes by adding the distance modulus for each star’s simulated heliocentric distance and the corresponding band extinction to its absolute magnitudes.

To reproduce the measurement uncertainties in each cluster, we derive empirical magnitude–uncertainty relations from all of its catalogue members (Hunt & Reffert 2024), rather than adopting general Gaia DR3 astrometric uncertainty models (Lindegren et al. 2021). Using locally weighted scatterplot smoothing (LOWESS), we fit the logarithmic uncertainties of the parallax and proper motions as functions of $G$, and the logarithmic uncertainty in each photometric band as a function of the magnitude in that band. For each mock star, we interpolate these relations at its noiseless magnitudes and add independent zero-mean Gaussian noise with the resulting standard deviations to its parallax, proper motions, and three magnitudes. We then retain sources with $4 \leq G \leq19$ mag.

We further restrict the sources to a footprint in sky position, parallax, and proper motion, defined from observed members of each cluster and fixed throughout the trials. The footprint radius is defined as $R_f = \max(3r_{95}, r_{95} + 20\ \mathrm{pc})$, where $r_{95}$ encloses 95 per cent of the observed members in projection at the catalogue distance $d_0$ (_dist50). The corresponding angular radius $\theta_f = \arctan(R_f/d_0)$ is limited to $0.5^\circ$–$10^\circ$. The proper motion cut is an ellipse centred on the catalogue proper motion, with semi-axes $2\Delta\mu$ in each component, where $\Delta\mu= \max(1.5\mu_{95}, \mu_{10}, 3e_\mu)$, limited to $0.3$–$5$ mas yr$^{-1}$, and $\mu_{10}$ corresponds to a transverse speed of $10\ \mathrm{km\ s}^{-1}$ at $d_0$. The parallax window is centred on the cluster parallax $\varpi_0$, with half-width $\Delta\varpi= \max(1.5\varpi_{95}, 3e_\varpi, \varpi_0R_f/d_0, 0.05\varpi_0)$, limited to $0.05$–$1$ mas. Here $\mu_{95}$ and $\varpi_{95}$ are the 95th percentiles of the members’ absolute offsets from the respective window centres, and $e_\mu$ and $e_\varpi$ are their median reported uncertainties, the term $\varpi_0R_f/d_0$ allows for the cluster’s line-of-sight depth. We retain stars within the angular radius, parallax window, and proper-motion ellipse:

Here, $\theta$ is the angular separation from the catalogue centre, and $\mu_{\alpha*,0}$ and $\mu_{\delta,0}$ are the catalogue proper motions.

We apply the same $4 \leq G \leq19$ mag and footprint cuts to the observed members and to the simulated cluster stars combined with mock field stars (§ 3.4), so that both counts refer to the same selection region. The selected simulated sample is then passed to HDBSCAN for membership identification (§ 3.5).

Mock Field Stars

With the synthetic cluster observations prepared, we add stars from the Gaia Early DR3 Mock Stellar Catalogue (GeDR3mock)⁵ (Rybizki et al. 2020). This catalogue was generated with galaxia (Sharma et al. 2011), combining a Besançon-like Galactic model with 3D extinction maps and stellar evolutionary models. It models the selection effects of the Gaia survey and provides Gaia-like uncertainty estimates.

For each cluster, we download field stars from a region with the same centre as the footprint and twice its angular radius $\theta_{\rm f}$, proper-motion scales $\Delta\mu$, and parallax half-width $\Delta\varpi$. Besides, we also apply additional quality criteria as described in Appendix B, those criteria are used only for download the field population, whereas the $4 \leq G \leq19$ mag and footprint cuts (§ 3.3) are applied afterwards to both populations. We retain the same field sample for each cluster throughout the calibration trials, providing a common background for the recovery comparison.

Membership Recovery and Initial-Mass Calibration

Membership recovery with HDBSCAN

We adopt HDBSCAN (Campello et al. [8]) because it was used to construct the homogeneous observational catalogue of Hunt & Reffert (2023, 2024) and is efficient for repeated forward calculations. The inputs are $(\mu_{\alpha*},\mu_{\delta},l,b,\varpi)$, with each dimension standardized by its mean and standard deviation. The main hyperparameters are the minimum cluster size $m_{\rm ClSize}$ and the number of samples in a neighbourhood $m_{\rm Pts}$. We fix $m_{\rm ClSize}=10$ and $m_{\rm Pts}=5$ and use excess of mass (eom) cluster selection, keeping recovery consistent across trials and limiting repeated tuning; the rationale is given in Appendix C1. These settings differ from those used in the all-sky search of Hunt & Reffert (2023) because our footprint selection and subsequent data standardization also differ. Following their procedure more strictly would allow a more rigorous comparison between the simulated and observed member samples. Our aim, however, is to develop an IMF-aware forward-modelling framework and derive IMF-aware selection functions, rather than reproduce the all-sky search.

From the HDBSCAN results obtained with the fixed settings, we identify the simulated target cluster as the non-noise label whose centroid is nearest to the observed cluster centre. Each source in this label receives a membership strength between 0 and 1, reported as a "probability". It is important to note that these values measure relative assignment strength rather than statistical membership probability and cannot be used to weight star counts; retaining all sources in the label would introduce disproportionately high field contamination. A suitable probability threshold must therefore be chosen: only sources that pass it are treated as simulated members in the final count comparison. The best threshold $P_{\rm best}$ is determined by the maximum difference between membership completeness and contamination:

where selecting sources in the target label with membership strength of at least $p$ gives the completeness $f_{\rm comp}(p)$ and contamination $f_{\rm cont}(p)$ (equation (6)). $P_{\rm best}$ is determined from simulated truth. Figure 2(a) illustrates this threshold selection for ASCC 105. This threshold is crucial because it controls the trade-off between completeness and contamination. It also changes the selected member count, which affects the convergence of mass calibration.

Vector (PDF)

At $P_{\rm best}$, we consider a simulated cluster successfully recovered if

where $N_{\rm sel}^{\rm sim}$ counts all selected sources, both correctly recovered members and field contaminants, while $N_{\rm cl}^{\rm sim}$ counts all cluster-origin particles passing the photometric and footprint cuts. We test convergence of the calibration only for recovered trials.

Initial cluster mass calibration

If the simulated cluster is recovered, we calibrate the trial against $N_{\rm sel}^{\rm obs}$, the observed member count after the same magnitude and footprint cuts, accepting differences of at most 20 per cent of this count. As illustrated in Figure 2(b) for ASCC 105, the selected count at $P_{\rm best}$ falls within this acceptance interval. To find a mass that meets this criterion, we start from $M_{\rm cl,0}=0.75\times N_{\rm sel}^{\rm obs}\ M_{\odot}$, following Zheng et al. (2026) and Marks & Kroupa (2012), and at each trial mass repeat the forward calculation, recovery, and threshold selection. We first search for a mass interval that brackets the observed count by doubling the trial mass when too few stars are selected, or halving it when too many are selected. We then bisect the bracketed mass interval; this assumes that $N_{\rm sel}^{\rm sim}$ generally increases with initial mass. Convergence is typically reached after a few simulations, although changes in the HDBSCAN recovery between trials can break this trend and prevent convergence.

If recovery fails, we instead compare counts after only the magnitude cut $4\leq G\leq19$: $N_{\rm mag}^{\rm sim}$ for simulated cluster stars with synthetic photometry, and $N_{\rm mag}^{\rm obs}$ for observed members. If $N_{\rm mag}^{\rm sim}<N_{\rm mag}^{\rm obs}$, the trial mass is too low and sets a lower bound for the bracketing interval, so we continue at higher masses; otherwise, we stop because the search direction is ambiguous. Such trials guide the search but are never accepted as calibrated models.

If mass calibration fails to converge with the fixed HDBSCAN settings, we repeat it with a search over 28 hyperparameter combinations at each trial mass: $m_{\rm ClSize}\in\{10,15,20,30,40,50,80\}$, $m_{\rm Pts}\in\{5,10\}$, and eom or leaf selection. In practice, only a small fraction of our sample required this additional search. For each combination, the target is the label with the most true cluster stars. We retain the combination maximizing the area under the completeness–contamination curve, then reselect $P_{\rm best}$ using equation (4). The search thus uses simulated truth for both target and hyperparameter selection; details are given in Appendix C2.

Membership Selection Functions

The selection functions compare the simulated samples before and after HDBSCAN membership recovery (§ 3.5.1). Before recovery, $N_{\rm all}^{\rm sim}$ counts all cluster-origin particles in the final snapshot (§ 3.2), including remnants and escapers; $N_{\rm all}^{\rm sim}$ counts those eligible for synthetic photometry and within the adopted magnitude limits; and $N_{\rm cl}^{\rm sim}$ counts the subset also inside the footprint, i.e. the cluster stars in the HDBSCAN input. At $P_{\rm best}$, the selected sample comprises correctly recovered cluster members and field contaminants, with counts $N_{\rm rec}^{\rm sim}$ and $N_{\rm cont}^{\rm sim}$, respectively, giving $N_{\rm sel}^{\rm sim}=N_{\rm rec}^{\rm sim}+N_{\rm cont}^{\rm sim}$. We define the membership selection functions as the magnitudedependent completeness $f_{\mathrm{comp}}(G)$ and contamination $f_{\mathrm{cont}}(G)$:

We use equal-frequency magnitude bins defined from the observed sample to accommodate their non-uniform magnitude distribution. Following the Rice rule, we set the number of bins to $k=\lceil2(N_{\mathrm{sel}}^{\mathrm{obs}})^{1/3}\rceil$. The global completeness and contamination are computed using the same ratios, with counts summed over the full adopted magnitude interval.

To improve the statistical robustness of our estimates, we perform 15 additional simulations for each accepted model, using different random seeds to generate the initial stellar masses, positions, and velocities while keeping the calibrated total mass and associated radius, IMF, age, observational-error model and noise seed, field population, and footprint fixed. Each realization is evolved independently and subjected to membership recovery, with $P_{\mathrm{best}}$ reselected but no mass recalibration. For models calibrated through a hyperparameter search, these repeats use the reference model’s final hyperparameters and target identification rule without further searches.

From the accepted model and 15 additional realizations, we retain only those meeting the global recovery criteria (equation (5)) and report the median and 16th–84th percentile range in each magnitude bin, along with the number of valid realizations. These summaries describe scatter among successfully recovered realizations under fixed model and field assumptions, rather than uncertainty in the median or systematic uncertainties.

Results

Magnitude-dependent completeness and contamination

We attempted to calibrate the initial mass of each of the 35 clusters for 3 assumed IMF slopes, $\alpha_2=-2.7,-2.3$, and $-1.9$, giving 105 combinations. The membership recovery procedure was run for 97 combinations: 92 with the fixed HDBSCAN hyperparameters and 5 through the supplementary hyperparameter search. The remaining 8 combinations failed the membership recovery criteria even after this search. The accepted models used in the following selection function analysis therefore cover 34 clusters, all except IC 2391, for which none of the three combinations satisfied the recovery criteria. Results are available for all three IMF slopes for 30 clusters, two slopes for three clusters, and one slope for IC 2602.

Across the full sample shown in Figure 3 and Fig. D1, median completeness typically lies between 0.8 and 1 at $G \lesssim17$, while median contamination is generally below 0.1. Towards the faint limit, completeness declines to $\sim0.5$–$0.8$ in many clusters, although others retain higher values. Contamination varies substantially between clusters, remaining below 0.1 in many faint bins but exceeding 0.5 in the faintest bins of several clusters, such as NGC 2232. These magnitude-dependent selection effects make corrections for incompleteness and contamination essential when deriving mass functions.

Vector (PDF)

Since the simulated clusters used to estimate these corrections depend on the assumed IMF, we examine whether and how the selection functions differ among the three values of $\alpha_2$. For about half of the clusters, the median completeness decreases steadily from $\alpha_2=-2.7$ to $-1.9$ over most of the magnitude range, as illustrated by ASCC 58, Collinder 140, and NGC 2232 in the left-hand column of Figure 3. However, the separations between the medians, typically $\lesssim0.1$, are comparable to or smaller than the 16th–84th percentile ranges across realizations, so the realization scatter does not allow us to establish whether this ordering is statistically significant. Moreover, even in these clusters, the shapes of the selection functions do not always follow the same trend for different $\alpha_2$: at $G\gtrsim17$, the curves for the three slopes cross, and the decline in completeness towards the faint limit differs between slopes. For example, in the faintest bin of Collinder 140, the completeness falls to 0.5 for $\alpha_2=-1.9$ but remains at 0.70–0.75 for the other two slopes. For the remaining clusters, the selection functions show no systematic dependence on $\alpha_2$, or the scatter is too large to identify one: the medians for the three slopes either nearly coincide, as for BH 99, or cross without a consistent order, as for Roslund 5 and IC 4665 in the right-hand column of Figure 3.

To assess whether these differences matter once the corrections are applied, we test how well they recover the true present-day luminosity function (PDLF) when derived with the correct IMF and with a different one. Because the selection functions are defined in magnitude bins, we perform this test on the luminosity function (LF) itself, without a mass–luminosity conversion. We treat each realization of the models with $\alpha_2=-1.9$ and $-2.7$ in turn as a target and correct its selected sample using the selection functions of other realizations of the same cluster. We compare three cases: no correction; same-IMF corrections, which use the other realizations of the target model; and wrong-IMF corrections, which use the realizations of the $\alpha_2=-2.3$ model. We use only realizations with valid counts in all magnitude bins, so that every LF is compared with the PDLF over the same full magnitude range.

In each magnitude bin $b$, let $N_b^{\mathrm{true}}$ be the true number of cluster stars in the target and $N_b$ the recovered number. Without correction, $N_b$ is the selected count, $N_b^{\mathrm{sel}}$; with correction, $N_b=N_b^{\mathrm{sel}}(1-f_{\mathrm{cont}})/f_{\mathrm{comp}}$, with $f_{\mathrm{comp}}$ and $f_{\mathrm{cont}}$ taken from one correction realization. We quantify the recovery with two diagnostics, evaluated for each pair of target and correction realizations and summarized for each cluster by the median over correction realizations and then over targets. The normalization bias,

is the signed relative error in the total number of stars, with positive values indicating overestimated counts. The shape distortion $L$ com… pares the LFs after normalizing each to unit total, so it is insensitive to the total number of stars,

where the median is taken in each bin before the absolute value. $L$ therefore measures the systematic change in LF shape, rather than the random scatter between realizations. For each cluster, we obtain 95 per cent intervals on $\epsilon_N$ and $L$ from 5000 bootstrap resamples of its target realizations, with the corrections kept fixed.

Figure 4 shows that the uncorrected LFs are biased in both normalization and shape. Without correction, $\epsilon_N$ is negative for almost all clusters, with medians of $-8$ and $-4$ per cent for targets with $\alpha_2=-1.9$ and $-2.7$, respectively, and the median $L$ is 2 and 1.5 per cent. Same-IMF corrections reduce both biases: the median $\epsilon_N$ becomes close to zero, the median $L$ decreases to 0.5 and 0.3 per cent, and the scatter between clusters decreases. Wrong-IMF corrections, which use the $\alpha_2=-2.3$ models, leave median $\epsilon_N$ values of $-1.4$ per cent for targets with $\alpha_2=-1.9$ and $+1.2$ per cent for those with $\alpha_2=-2.7$. Most of the normalization bias is therefore still removed, but the total counts remain underestimated when the target IMF is flatter than $\alpha_2=-2.3$ and overestimated when it is steeper. The shape distortion, however, is reduced much less. Figure 4(c) and (d) show that wrong-IMF corrections leave median $L$ values of 1.6 and 1.1 per cent, about three times the same-IMF values and only modestly below the uncorrected ones. Besides this higher median, $L$ also varies more between clusters than with same-IMF corrections, suggesting that clusters respond differently to the IMF mismatch.

Vector (PDF)

Overall, corrections derived with the correct IMF largely recover both the normalization and the shape of the LF, whereas those derived with a mismatched IMF bring the total number of stars close to its true value but still leave a shape distortion.

Global completeness and contamination

Although the selection functions in Figure 3 do not show a clearly resolved ordering with $\alpha_2$ in individual magnitude bins, such an ordering becomes systematic once the counts are integrated over the full magnitude range. Figure 5 shows the resulting global completeness and contamination for the 30 clusters with accepted models at all three IMF slopes. Each point in Figure 5(a) is the median over valid realizations for one cluster at one slope. At fixed $\alpha_2$, the clusters span a wide range, from 0.77 to 0.98 in completeness, with interquartile ranges of 0.05–0.06. Most grey lines linking the same cluster nevertheless decline in completeness and rise in contamination from $\alpha_2=-2.7$ to $-1.9$.

Vector (PDF)

To separate the effect of the IMF from these differences between clusters, we compute, for each cluster, the difference between the values at two slopes, defined as the value at the flatter slope minus that at the steeper one. Figure 5(b) shows these differences and their median across clusters. We estimate a 95 per cent interval on each median from 20 000 bootstrap resamples of the clusters, using the bias-corrected and accelerated (BCa; (Efron 1987)) method. Between $\alpha_2=-2.7$ and $-1.9$, the median difference is $-0.040$ in completeness, with a BCa interval of $-0.053$ to $-0.031$, and $+0.015$ in contamination, with an interval of 0.008 to 0.021. The completeness is lower at $\alpha_2 = -1.9$ in all 30 clusters, and the contamination is higher in 29 of them. The two intermediate steps show smaller shifts in the same direction, $-0.012$ and $-0.024$ in completeness, with decreases in 21 and 28 clusters, respectively. Thus, even though the models for each cluster are calibrated against the same observed member count, those with flatter IMFs systematically yield lower global completeness and higher global contamination. When the same observed members are corrected for these selection effects, a flatter IMF therefore implies a larger upward correction of the star counts for incompleteness and a larger downward correction for contamination. This is consistent with the signs of the residual normalization bias for wrong-IMF corrections in Figure 4(a) and (b).

Although systematic, this IMF-dependent offset is modest compared with the differences between clusters at fixed $\alpha_2$. The median shift of 0.04 in completeness between the two extreme slopes is comparable to the interquartile range across clusters but much smaller than the full range, and the shift of 0.015 in contamination is likewise comparable to its interquartile range of 0.01–0.02. Moreover, the clusters largely retain their relative order across the three slopes, as indicated by the predominantly parallel grey lines in Figure 5(a). We next examine how the differences between clusters relate to their properties.

Figure 6(a) shows that older simulated clusters are recovered less completely and with more contamination. The Spearman coefficients

Vector (PDF)

$(\rho)$ with age are $-0.87$, $-0.85$, and $-0.82$ for completeness and 0.60, 0.50, and 0.70 for contamination, for $\alpha_2 = -2.7$, $-2.3$, and $-1.9$, respectively; all six 95 per cent BCa intervals exclude zero. However, age only sets how long each model has evolved. Clusters of the same age can be at very different stages of dynamical evolution, depending on their mass and size, so age alone does not identify the underlying physical mechanism.

We therefore further consider a dynamical clock, the cumulative number of half-mass relaxation times, $N_{\mathrm{relax}} = \int_0^t \mathrm{d}t'/t_{\mathrm{rh}}(t')$, where $t$ is the cluster age and $t_{\mathrm{rh}}(t')$ is the half-mass relaxation time (Spitzer 1987) of the simulated cluster at time $t'$. For $\alpha_2 = -2.7$ and $-2.3$, completeness decreases with $N_{\mathrm{relax}}$ about as strongly as with age, with Spearman coefficients of $-0.82$ and $-0.88$. This is consistent with cumulative two-body relaxation contributing to the structural evolution that makes the simulated clusters progressively more difficult to recover. For $\alpha_2 = -1.9$, however, the coefficient is only $-0.29$, with an interval of $-0.60$ to 0.10, and the contamination coefficients are also weaker. One possible reason for the weaker relation⁶ at $\alpha_{2} = -1.9$ is that $N_{\mathrm{relax}}$ measures the progress of two-body relaxation, whereas a flatter IMF contains more massive stars, whose stellar evolution mass loss may also contribute to the early expansion of the cluster. This expansion is not driven by two-body relaxation; moreover, the increase in $r_{\mathrm{h}}$ lengthens $t_{\mathrm{rh}}$, reducing the number of relaxation times accumulated by the same age. For the flatter IMF, $N_{\mathrm{relax}}$ may therefore be a less reliable proxy for the cluster structure relevant to membership recovery. This suggests that $N_{\mathrm{relax}}$ alone may not consistently explain differences in recovery between clusters across the three IMF assumptions.

Rather than relying on age or a dynamical clock as proxies, we next examine the structural evolution of the simulated clusters directly. Figure 6(c) uses the ratio $r_{\mathrm{h}}/r_{\mathrm{h},0}$ of the half-mass radius at the cluster age to its initial value. More strongly expanded models show lower completeness for all three IMF slopes, with coefficients of $-0.89$, $-0.88$, and $-0.85$, and higher contamination, with coefficients of $0.45$–$0.78$. Unlike $N_{\mathrm{relax}}$, the expansion remains closely related to membership recovery even for $\alpha_{2} = -1.9$. This suggests that the structural state of the model, as traced by its expansion, is more consistently associated with recovery than the dynamical clock alone.

Finally, Figure 6(d) uses the projected stellar surface density $\Sigma_{50} = 0.5N_{\mathrm{cl}}^{\mathrm{sim}}/(\pi R_{50}^{2})$, where $R_{50}$ is the projected radius enclosing half of the simulated cluster members. Whereas $r_{\mathrm{h}}/r_{\mathrm{h},0}$ shows how much a model has expanded relative to its initial state, $\Sigma_{50}$ measures how compact it is at the cluster age. Completeness increases with $\Sigma_{50}$, with coefficients of $0.93$, $0.91$, and $0.84$, and contamination decreases, with coefficients of $-0.68$, $-0.49$, and $-0.72$. The absolute coefficients are comparable to or slightly larger than those for age and expansion, although their confidence intervals overlap. These trends are consistent with the density-based nature of HDBSCAN: a more diffuse stellar distribution provides less density contrast against the field, making cluster members more difficult to separate from surrounding sources. $\Sigma_{50}$ may therefore serve as a link between the evolution traced in Figure 6(a)–(Figure 6c) and the recovery performance.

Taken together, the four panels of Figure 6 show that recovery is strongly related to age, expansion, and surface density for all three IMF slopes, but not uniformly to $N_{\mathrm{relax}}$. These relations are consistent with a picture in which models expand and become less dense as they evolve, making their members harder to separate from field stars. Because the four quantities are correlated with one another, however, the present analysis cannot isolate the effect of any single one.

Discussion

Effect of the assumed IMF

We now examine whether the systematic completeness offset between the different IMF assumptions is accompanied by changes in the evolution and structure of the simulated clusters. For each cluster, Figure 7 shows the changes from $\alpha_2=-2.7$ to $-1.9$ in stellar evolution mass loss, structure, and global completeness. Because both models of a cluster share its age, orbit, field population, and observed member count, these changes remove the offsets between clusters seen in Figure 6, although they still include the effect of the different initial masses calibrated for the two IMF slopes.

Vector (PDF)

We quantify the stellar evolution mass loss by $f_{\mathrm{SE,loss}}=1-M_{\mathrm{all}}(t)/M_{\mathrm{cl},0}$, the cumulative fraction of the initial cluster mass lost by the cluster age $t$ through stellar and binary evolution, including mass ejected in mergers; $M_{\mathrm{all}}(t)$ is the total mass of all particles in the simulation at $t$, including remnants and escapers, so that stars escaping from the cluster are not counted as mass loss. As shown in Figure 7(a), in 27 of the 30 clusters the flatter-IMF model has a larger $f_{\mathrm{SE,loss}}$, with $\Delta f_{\mathrm{SE,loss}}\approx0.15$–$0.4$, and expands more relative to its initial half-mass radius, with $\Delta\log(r_{\mathrm{h}}/r_{\mathrm{h},0})\approx0.25$–$0.65$. A flatter IMF places a larger fraction of the cluster mass in massive stars, leading to stronger stellar evolution mass loss and, consequently, greater cluster expansion (e.g. [36] (Haghi et al. 2020)). However, the amount of additional mass loss does not scale monotonically with the amount of additional expansion, because the resulting dynamical response also depends on how rapidly and where the mass is lost (e.g. [1] (Hills 1980)).

Regardless of the detailed origin of the expansion, its effect on the present-day structure is clear. Clusters for which the flatter-IMF model expands more strongly also show a systematically larger decrease in $\Sigma_{50}$, as in Figure 7(b). The Spearman coefficient is $-0.90^{+0.18}_{-0.06}$ and most clusters lie along a narrow sequence. Thus, the differences in expansion between the two IMF assumptions are reflected in the projected stellar densities of the stars passed to membership recovery.

Figure 7(c) shows that larger decreases in $\Sigma_{50}$ tend to accompany larger decreases in completeness, and Figure 7(d) shows a similar tendency for models that expand more. However, the evidence for these trends is limited. The Spearman coefficients are 0.37 and $-0.41$, with 95 per cent intervals of $-0.04$ to 0.66 and $-0.71$ to 0.00, respectively. For $\Sigma_{50}$, this is likely because HDBSCAN recovery depends on the density contrast between cluster and field stars in the five-dimensional space $(l,b,\varpi,\mu_{\alpha*},\mu_\delta)$, which a projected cluster density alone cannot describe.

Three clusters, ASCC 19, OC 0280, and Alessi–Teutsch 10, have $\Delta f_{\mathrm{SE,loss}}$, $\Delta\log(r_{\mathrm{h}}/r_{\mathrm{h},0})$, and $\Delta f_{\mathrm{comp}}$ all close to zero in Figure 7. They are younger than 8 Myr and have the lowest calibrated initial masses in our sample, below 60 $M_\odot$ under both IMF assumptions. Under optimal sampling, these low initial masses imply low maximum stellar masses through equation (2). For the resulting stellar mass distributions, little stellar evolution mass loss has occurred by the cluster ages, even with the flatter IMF. By contrast, other clusters of similar age but higher initial mass already show $\Delta f_{\mathrm{SE,loss}}\approx0.15$–$0.18$ and $\Delta\log(r_{\mathrm{h}}/r_{\mathrm{h},0})\approx0.5$. This contrast suggests that changes in mass loss and expansion between the two IMF assumptions depend on cluster mass as well as age.

Overall, Figure 7 shows that models with $\alpha_2=-1.9$ generally undergo more stellar evolution mass loss and expand more than those with $\alpha_2=-2.7$. Their stronger expansion is closely reflected in lower projected stellar densities, so these models are more diffuse. The systematically lower global completeness at $\alpha_2=-1.9$ is then consistent with more diffuse clusters being harder to recover. However, clusters with similar structural changes can still differ considerably in their completeness shifts, suggesting that structure is not the only factor behind the completeness difference between the two slopes.

Algorithm dependence and applicability of the framework

Our selection functions are specific to HDBSCAN with the settings described in Section 3.5. Open cluster member lists, however, are constructed with a variety of methods, including DBSCAN (Castro-Ginard et al. 2020, 2022), fastMP (Perren et al. 2023), pyUPMASK (Pera et al. 2021; Qin et al. 2023), STARGO (Yuan et al. 2018; Pang et al. 2022), and other machine-learning approaches (Van Groeningen et al. 2023). Because these methods identify members by different criteria, the resulting lists would differ in their completeness at the faint end and in their levels of field contamination. Even with HDBSCAN, changing the hyperparameters or probability threshold can alter the recovered sample. Selection functions must therefore be derived for the specific procedure used to construct each member list.

The forward-modelling framework developed in this work extends naturally to these algorithm-specific selection functions. Of the four steps in Figure 1, only the membership recovery in step 3 depends on the algorithm. To derive the selection functions of a given member list, step 3 should reproduce the procedure used to construct that list, applied to the simulated sample with the same input variables, preprocessing, hyperparameters, and member selection criteria as for the observed data. The observational inputs to steps 1 and 2 are the initial cluster mass $M_{\mathrm{cl},0}$ or are then taken from, or recalibrated against, that list. Thus, although the selection functions are specific to one algorithm, the framework from which they are derived can be generalized to others.

Because the origin of each simulated star is known, the framework can also be used to evaluate how well candidate algorithms or settings recover a given cluster, and thus to choose the procedure for determining its membership. Given the close relation between HDBSCAN recovery and cluster structure in § 4.2, such evaluations should use simulations of the target cluster rather than a generic mock sample.

Caveats

The selection functions derived here depend on several model choices. We vary only the IMF slope above 0.5 $M_\odot$, keeping the slope below this mass fixed and assigning stellar masses by optimal sampling, so that the differences between the models of a cluster can be attributed to $\alpha_2$. The slope below 0.5 $M_\odot$, which governs the number of the faintest members, is therefore not tested. Optimal sampling also fixes the maximum stellar mass for a given $M_{\mathrm{cl},0}$ through equation (2). In the three clusters with the lowest $M_{\mathrm{cl},0}$, this mass is too low for the models with $\alpha_2=-2.7$ and $-1.9$ to differ appreciably in stellar evolution mass loss, expansion, or completeness, as shown in § 5.1. Stochastic sampling could add more massive stars to such clusters, making their selection functions more sensitive to $\alpha_2$, and would increase the scatter between realizations.

The initial half mass radii follow the relation of Marks & Kroupa (2012), which describes embedded clusters before gas expulsion, whereas our models start gas-free and in virial equilibrium, without primordial binaries or mass segregation. They may therefore remain more compact than the observed clusters (Zheng et al. 2026), many of which are expanding (Della Croce et al. 2024). As completeness increases with $\Sigma_{50}$ in Figure 6(d), such models would tend to overestimate the completeness, and the offset between $\alpha_2=-2.7$ and $-1.9$ found in § 5.1 may also depend on the adopted initial structure. Gas expulsion and primordial binaries can be added to the $N$-body models without changing the rest of the framework.

As noted in § 3.5.1, our membership recovery follows the

search used to construct the observed member lists only approximately, and $P_{\mathrm{best}}$ is chosen using simulated truth, which is unavailable for observed data. The selection functions therefore describe recovery with an optimally chosen threshold rather than with the criteria of the observed lists. Reproducing the member search of the observed catalogue in full, with the same input variables, preprocessing, hyperparameters, and member selection criteria, would give selection functions that apply directly to the observed lists.

Finally, the calibration matches only $N_{\mathrm{obs}}^{\mathrm{sel}}$. This fixes $M_{\mathrm{cl},0}$ for each assumed IMF but cannot distinguish between IMFs or uniquely determine the initial structure. The choice of IMF nevertheless matters: Figure 4(c) and (d) show that corrections derived with $\alpha_2=-2.3$ leave residual distortions in the LF shapes of the target models with $\alpha_2=-1.9$ and $-2.7$. When correcting an observed cluster, the selection functions should therefore be inferred jointly with its IMF, rather than taken from a model with an IMF assumed in advance. Such joint inference requires observables beyond the member count. Matching the luminosity function and density profile of the selected members would constrain $\alpha_{2}$ and $r_{\mathrm{h},0}$, respectively; with a larger grid of $N$-body models, the framework could then infer these initial conditions together with $M_{\mathrm{cl},0}$ and the corresponding selection functions.

Conclusions

In this work, we developed an IMF-aware forward-modelling framework to quantify the membership selection functions—namely, completeness and field star contamination—of open clusters. We evolve cluster models with $N$-body simulations (Wang et al. 2020b), convert their outputs into synthetic Gaia observations, and add mock field stars (Rybizki et al. 2020). We then recover cluster members with HDBSCAN (Campello et al. [8]; Hunt & Reffert 2023) to quantify membership completeness and contamination, examining how these quantities vary with the assumed stellar initial mass function (IMF) and relate to cluster evolution.

Applying our framework to 35 nearby young open clusters with three assumed high-mass ($\ge0.5M_{\odot}$) IMF slopes ($\alpha_{2} = -2.7$, $-2.3$, $-1.9$), we find the following main results:

• Magnitude- and cluster-dependent selection: Membership completeness and contamination are highly magnitude-dependent and vary between individual clusters. Faint members are systematically harder to separate from field stars, confirming that accurate derivations of the present-day mass function require corrections that are specific to both the cluster and the precise membership identification procedure used.

• IMF dependence of selection functions: The global completeness and contamination depend systematically on the assumed IMF. Models initialized with a flatter IMF ($\alpha_{2} = -1.9$) contain a larger fraction of massive stars, leading to stronger stellar evolution mass loss. Consequently, these clusters expand more ([36] (Chernoff & Weinberg 1990; Haghi et al. 2020; Wang & Jerabkova 2021)) and exhibit lower projected stellar surface densities at the present day. This reduced density contrast against the field makes members more difficult to recover, yielding systematically lower global completeness and higher global contamination compared to models with steeper IMFs.

• Impact on the recovered luminosity function: Applying completeness and contamination corrections derived from a mismatched IMF can largely remove the normalization bias (recovering the total number of stars), but it leaves residual distortions in the shape of the present-day luminosity function.

• Correlation with structural evolution: Cluster recovery performance is strongly tied to the structural state of the model. Completeness significantly decreases and contamination increases for models that have expanded more relative to their initial radii and for those with lower present-day projected stellar surface densities, consistent with the density-based nature of HDBSCAN.

Overall, our results suggest that membership selection functions cannot be universally assumed or derived from generic mock clusters. Because a cluster’s evolution and present-day structure are sensitive to its initial mass distribution, correcting an observed cluster’s mass function using selection functions derived from a mismatched IMF would inherently bias the result. To accurately infer the IMF of an open cluster, these membership selection effects should ultimately be inferred jointly with the IMF and other initial cluster conditions, by incorporating membership selection into $N$-body inferences of these conditions ((Baumgardt et al. 2023; Weis et al. 2025; Hobart et al. 2026; Zheng et al. 2026)).

Acknowledgements

The authors appreciate the valuable comments and suggestions from the reviewer. This work is supported by the National Key Research and Development Program of China No. 2025YFF0511001. L.W. thanks the support from the National Natural Science Foundation of China through grant 12573041 and 12233013, the GuangDong Basic and Applied Basic Research Foundation (2026A1515012460), the High-level Youth Talent Project (Provincial Financial Allocation) through the grant 2023YHSPT0706, the Fundamental Research Funds for the Central Universities, Sun Yat-sen University (2025QNPY04). X.P. acknowledges the financial support of the National Natural Science Foundation of China through grants 12573036 and 12233013, and the China Manned Space Program with grant No. CMS-CSST-2025-A08. H.T. is supported by the National Key R&D Program of China No. 2025YFF0511002 and 2024YFA1611902 and the China Manned Space Project with grant No. CMS-CSST-2025-A11. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Data availability

Data are available upon request.

Cluster identifications and adopted properties

Table A1 lists the identifications and properties of the 35 clusters in the observational sample described in § 2. All numerical values are taken from the catalogue of Hunt & Reffert (2024); the names used by Pang et al. (2022) are included to document the sample selection. We match these names against the AllNames column of the Hunt & Reffert (2024) catalogue. Where multiple candidate entries are found, we adopt the entry with the smallest angular separation, absolute distance difference, and proper motion difference from the corresponding object in Pang et al. (2022).

Ages from Hunt & Reffert (2024) are expressed in linear units by converting the catalogue columns logAge16, logAge50, and logAge84 using $t_q/\mathrm{Myr}=10^{\ell_q-6}$, where $\ell_q$ denotes the corresponding logarithmic age and $q \in\{16,50,84\}$. They are reported as t_{50}^{+(t_{84}-t_{50})}_{-(t_{50}-t_{16})}. The simulations use $t_{50}$ rounded down to the nearest Myr.

The member count $N_{\mathrm{cat}}^{\mathrm{obs}}$ is taken from column N of the Hunt & Reffert (2024) catalogue, before the magnitude and observational footprint cuts described in § 2. Cluster masses are taken directly from the MassTot column, whose values already include corrections for completeness and unresolved binaries. The quoted mass uncertainties are taken from e_MassTot. These are present-day total masses, rather than initial cluster masses or masses enclosed within the Jacobi radius.

GeDR3mock query for field stars

Table B1 lists the conditions used to retrieve the mock field stars of each cluster from gedr3mock::main through the TAP service of the GAVO Data Center (https://dc.g-vo.org/tap).

HR24: Hunt & Reffert (2024); P22: Pang et al. (2022). Masses and their uncertainties are taken directly from MassT oT and e_MassT oT, respectively; the masses already include corrections for completeness and unresolved binaries. $N_{\mathrm{cat}}^{\mathrm{obs}}$ is the HR24 catalogue member count (column N) before the magnitude and observational footprint cuts applied in this work. Superscript m denotes an HR24 moving-group classification.

HDBSCAN hyperparameter settings

Fixed hyperparameters

By default, we adopt $m_{\mathrm{clSize}} = 10$, $m_{\mathrm{Pts}} = 5$, and excess of mass (eom) cluster selection for all models. Holding these settings fixed keeps the recovery procedure consistent across mass-calibration trials and limits repeated tuning.

The neighbourhood parameter $m_{\mathrm{Pts}}$ sets the scale for local density estimation. Larger values make clustering more conservative, generally assigning more sources to noise and restricting detections to denser regions⁷. We choose $m_{\mathrm{Pts}} = 5$ as a relatively permissive setting; this value has also been used for open cluster membership identification by Alfonso et al. (2024). The minimum cluster size $m_{\mathrm{clSize}} = 10$ allows groups with as few as ten sources to remain eligible for recovery, as also adopted by Qin et al. (2023).

The eom method selects structures by their integrated stability in the density hierarchy, while leaf selection produces finer subdivisions. Our aim is to recover the target cluster as a whole, motivating the use of eom. These choices define the recovery procedure for which the selection functions are measured. Each trial must satisfy the recovery criteria in § 3.5.1; if the mass calibration does not converge with these settings, we perform the supplementary hyperparameter search below.

Quantity	Condition
Search limits
Sky position	$\theta\le\theta_{\mathrm{d}}$
Parallax	$\lvert\varpi-\varpi_0\rvert\le\Delta\varpi_{\mathrm{d}}$
Proper motion	$\sum_i [(\mu_i-\mu_{i,0})/\Delta\mu_{i,\mathrm{d}}]^2 < 4$
Magnitude	phot_g_mean_mag $\le21$ mag
Quality criteria
Visibility periods	visibility_periods_used $> 5$
Population	popid $\ne11$
$G$ flux	phot_g_mean_flux_over_error $> 10$
$G_{\mathrm{BP}}$ flux	phot_bp_mean_flux_over_error $> 10$
$G_{\mathrm{RP}}$ flux	phot_rp_mean_flux_over_error $> 10$
Parallax	parallax_over_error $> 10$

Free hyperparameter search

The supplementary search covers 28 combinations of $m_{\mathrm{clSize}} \in\{10, 15, 20, 30, 40, 50, 80\}$, $m_{\mathrm{Pts}} \in\{5, 10\}$, and eom or leaf selection. For each fit, the target is the non-noise label containing the most true cluster stars; ties favour the smaller label.

We evaluate completeness and contamination at 500 equally spaced membership-strength thresholds over $[0,1]$, sort by contamination, and integrate completeness against contamination using the trapezoidal rule. Empty selections are assigned zero contamination for ranking only. We retain the fit with the largest area; ties favour eom over leaf, then increasing $m_{\mathrm{clSize}}$ and $m_{\mathrm{Pts}}$. This completeness–contamination area is not a standard receiver-operating-characteristic area. We then select $P_{\mathrm{best}}$ using the trial-wide $f_{\mathrm{comp}}-f_{\mathrm{cont}}$ criterion in § 3.5.1, without refitting. The supplementary search therefore uses simulated truth for both target and hyperparameter selection, independently of the observed member count.

Selection functions of the full sample

Fig. D1 shows the selection functions of the 28 clusters with accepted models that are not included in Figure 3; together, the two figures cover all 34 clusters used in the selection function analysis. Accepted models are available for only two IMF slopes for NGC 7429, NGC 2451A, and NGC 1980, and only for $\alpha_2=-2.7$ for IC 2602, for which no comparison between slopes is therefore possible. For IC 2602, only three of the 16 realizations satisfy the recovery criteria in equation (5), so its percentile ranges are poorly sampled.

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)