© The Authors 2024

1 Introduction

Stellar clusters (SCs) are the building blocks of our and other galaxies (Kroupa 1995a; Kroupa & Boily 2002; Lada & Lada 2003; Kroupa 2005; Portegies Zwart et al. 2010). Therefore, understanding their masses and the mass functions (MFs) within them is of the utmost importance for comprehending the Universe as we see it today. One factor that strongly influences the measurements of these properties are unresolved binaries (Kroupa et al. 1991; Kouwenhoven & de Grijs 2008; Gieles et al. 2010; Kroupa & Jerabkova 2018; Borodina et al. 2019; Rastello et al. 2020).

While the binary fraction in the Galaxy is ≈0.5 (Duquennoy & Mayor 1991; Raghavan et al. 2010; Marks & Kroupa 2011, and references therein), SCs are likely born with up to 100% of their stars in binaries (Mathieu 1994; Kroupa 1995a,b). The difference between the observed distribution functions of binary systems in the Galactic field in comparison to that observed in star-forming regions allowed Kroupa (1995a,b) to deduce the properties of the initial distribution functions of young binary-rich populations. There is still some debate around whether this is also true for very low-mass binaries, many of which might dissolve before reaching the main sequence or even not start out as binaries (Thies et al. 2015; Offner et al. 2023). Due to the dense environment of the clusters, the initially high binary fractions in SCs are reduced quickly (Kroupa 1995a,b; Ivanova et al. 2005; Goodwin & Kroupa 2005; Kobulnicky & Fryer 2007). The stars lost during SC evolution form the field stars observed in the Galaxy today (Kroupa 2002b; Marks & Kroupa 2011; Dabringhausen et al. 2022; Ferrone et al. 2023).

Dinnbier & Kroupa (2020) categorised these stars into tail 1 and tail 2 (T1 and T2) classes. T1 stars are lost from the birth-embedded cluster due to the decrease in the gravitational potential during gas expulsion. The T2 stars, on the other hand, evaporate slowly from the SC after this initial phase. Both regimes overlap in phase space and so are not observed as two separated regions. They are also contaminated by high-velocity stars ejected due to few-body interactions and supernova (SN) kicks. These tails are expected to inherit binaries from the SC, and if the binaries are unresolved, this has implications for the observations of the tails and the SC.

One of the most common ways to determine the mass of a star is to use the measured luminosity. Many studies have determined the mass-luminosity relation (MLR) for main-sequence stars in different mass ranges (see e.g. Hertzsprung 1923; Russell et al. 1923; Kroupa et al. 1993; Henry & McCarthy 1993; Eker et al. 2015, and references within). However, if a binary is unresolved, only the combined luminosity of its members can be measured. An unassuming observer would take this combined luminosity and use it to compute the mass of the assumed single star. This mass, however, does not equal the sum of the masses of the binary components. Therefore, an observer would arrive at a wrong mass for the binary and the entire SC.

Alternatively, the mass of a virialised SC can be deduced dynamically because its dynamical mass is proportional to the square of the line-of-sight velocity dispersion (e.g. Illingworth 1976; Spitzer 1987). However, for this deduction, a binary fraction of 0 is assumed. While this mass determination might still work if centre-of-mass velocities for binaries are used, many studies have shown that the velocity dispersions for systems with large fractions of unresolved binaries are overestimated (see e.g. Kouwenhoven & de Grijs 2008; Gieles et al. 2010; Borodina et al. 2019; Rastello et al. 2020; Borodina et al. 2021). Therefore, binaries will cause an overestimate of the SC mass if it is measured dynamically. However, long-duration observations over multiple epochs can mitigate the effect (Cottaar et al. 2012; Ramírez-Tannus et al. 2024). Cottaar & Hénault-Brunet (2014) suggest that two epochs are sufficient for this in many SCs; however, they also point out that many binaries have periods of thousands of years, making it practically impossible to detect them.

Similarly, the presence of binaries influences the deduced MF for low-mass stars (Kroupa et al. 1991). Having an unresolved binary means that, when determining the MF, two lower-mass stars are missing and one slightly more massive star is added (Kroupa & Jerabkova 2018). This makes the MF appear more top-heavy. Li et al. (2023) investigated this effect using mock data, and it was further investigated using N-body simulations in this work.

The purpose of this work is, therefore, to determine how the binary fraction affects our understanding of the mass and MF in SCs and their tails. For the first time, this is done using realistic N-body simulations of SCs to determine the development of the binary population. Under the assumption that all binaries are unresolved, possible conclusions an observer would draw from these data are computed and compared to the results obtained under the assumption that all binaries are resolved. This work is structured as follows: Section 2 explains the methods used for this work. The results are discussed in Sect. 3, and Sect. 4 concludes the paper.

2 Methods

2.1 Description of the simulations

The SC models initially follow non-mass-segregated Plummer density profiles and the Marks & Kroupa (2012) relation between the initial cluster mass in stars (M_ecl) and the half-mass radius (r_h). This relation was deduced from the density at birth of open clusters required to account for the observed binary populations in open clusters. It is important to emphasise that the SCs are non-mass-segregated, because the stars in the outer regions of the SCs are preferentially lost during gas expulsion and through evaporation. On the other hand, dynamical ejection is more important in the SC centre, as the high density of stars leads to a large number of encounters. Since evaporation and gas expulsion are the main ways stars are lost, a mass-segregated SC would lead to the tidal tails becoming bottom-heavy, while the SC would become bottom-light.

Table 1 lists the initial parameters of the main model, M1, used in this work and a number of lower-mass models, initialised with different random seeds, used for additional analysis. Initially, the clusters are embedded in their natal gas, which is approximated by an analytical Plummer model of the same halfmass radius as the stellar population. The total gaseous mass is 2 M_ecl in each model; this assumes a star formation efficiency of 1/3 (see e.g. Megeath et al. 2016).

Up to the age of 0.6 Myr, the gaseous component is constant in time; after that, its density decreases exponentially on a timescale of r_h/(10 km s⁻¹), which approximates feedback from young massive stars. More details of the gas expulsion model can be found in Kroupa et al. (2001, see also Lada et al. 1984 and Goodwin 1997 for earlier analytic models of gas expulsion and Calura et al. 2015 and Dinnbier & Walch 2020 for hydrodynamical simulations). The clusters move on circular orbits within the Galactic disc at galactocentric radius R_g = 8 kpc; the Galactic potential is represented by the Allen & Santillan (1991) model.

Initial positions and velocities of stars correspond to virial equilibrium with their primordial gas; they were generated using the Aarseth et al. (1974) method. Individual stellar masses were sampled from the two-part-power-law canonical initial mass function (IMF; Kroupa 2001; Kroupa et al. 2013): $$\xi (m)=\{\begin{array}{ll}{k}_1{m}^{-1.3}& m<0.5M_{\odot },\\ k_2{m}^{-2.3}& m\ge 0.5M_{\odot }.\end{array}$$(1)

Here k₁ and k₂ are constants chosen such that the IMF is continuous and the integration over the IMF yields the total mass of the SC. The stars have masses m > 0.08 M_⊙, brown dwarfs are neglected as they constitute a negligible mass contribution.

All stars are paired to binaries. To match observational quantities (e.g. Duquennoy & Mayor 1991; Sana et al. 2012; Kobulnicky et al. 2014), the distributions for mass ratio, orbital period, and eccentricity are separated into a set used for lowermass (primary mass m₁ < 3 M_⊙) and one for massive (m₁ ≥ 3 M_⊙) binaries (Kroupa 1995a,b; Belloni et al. 2017). We note that the boundary between high- and low-mass stars is at 5 M_⊙ in Belloni et al. (2017), while in the current work it is 3 M_⊙. This has a negligible effect on the circularisation and secondary masses of stars in this range. To compute the initial stellar properties, first the pre-main-sequence stars are generated following (Kroupa 1995b) with the updates in Belloni et al. (2017):

The masses of the primary, m₁, and the secondary, m₂, are drawn randomly from the Kroupa (2001) IMF. For the highmass stars (m₁≥3 M_⊙), stars are ordered before being paired: the most massive star is paired with the second most massive one, and so on. The low-mass stars (m₁ < 3 M_⊙) are paired randomly. The fraction of orbits with an eccentricity between e and de, dN = f_ede, is described by Eq. (4b) of Kroupa (1995a): $$f_e=2e,$$(2)

and the fraction of periods between log₁₀ P and log₁₀ P + d log₁₀ P, dN = f_Pd log₁₀ P is described by Eq. (8) of Kroupa (1995b): $$f_p=2.5\frac{{\log }_{10}(P/\text{ days })-1}{45+{[{\log }_{10}(P/\text{ days })-1]}^2}.$$(3)

For the massive stars, the pre-main-sequence phase is assumed to be too short to change the properties of the binaries significantly. However, for the low-mass stars, the properties are changed due to pre-main-sequence eigenevolution as computed by Kroupa (1995b) and summarised by Belloni et al. (2017): $$\rho ={\left(\frac{28R_{\odot }}{R_p}\right)}^{\frac{3}{4}},$$(4) $$\ln \left(e_{\mathrm{i}\mathrm{n}\mathrm{i}}\right)=-\rho +\ln \left(e_{\mathrm{b}\mathrm{i}\mathrm{r}}\right)$$(5) $$q_{\mathrm{i}\mathrm{n}\mathrm{i}}=\{\begin{array}{ll}{q}_{\mathrm{b}\mathrm{i}\mathrm{r}}+\rho (1-q_{\mathrm{b}\mathrm{i}\mathrm{r}}),& \rho \le 1,\\ 1,& \rho >1,\end{array}$$(6) $$m_{1,\mathrm{i}\mathrm{n}\mathrm{i}}=m_{1,\mathrm{b}\mathrm{i}\mathrm{r}},$$(7) $$m_{2,\text{ini}}=q_{\text{ini }i}{m}_{1,\text{bir}},$$(8) $$P_{\mathrm{i}\mathrm{n}\mathrm{i}}=P_{\mathrm{b}\mathrm{i}\mathrm{r}}\sqrt{\frac{m_{1,\mathrm{b}\mathrm{i}\mathrm{r}}+m_{2,\mathrm{b}\mathrm{i}\mathrm{r}}}{m_{1,\mathrm{i}\mathrm{n}\mathrm{i}}+m_{2,\mathrm{i}\mathrm{n}\mathrm{i}}}}{\left(\frac{1-e_{\mathrm{b}\mathrm{i}\mathrm{r}}}{1-e_{\mathrm{i}\mathrm{n}\mathrm{i}}}\right)}^{\frac{3}{2}}.$$(9)

Here the subscript ‘ini’ denotes that a variable is at the end of the eigenevolution during the pre-main-sequence time, while the subscript ‘bir’, for birth value, is used for variables at the beginning of the pre-main-sequence phase. R_p is the pericentre distance of the binary, q ≤ 1 the fraction between secondary mass and primary mass, m₁ and m₂ are the primary and secondary mass, respectively, and P the period of the binary. We note that the secondary masses are slightly increased through pre-main-sequence eigenevolution, flattening the IMF negligibly.

In this work, we define the binary fraction as $f_{\text{bin }}=\frac{N_{\text{bin }}}{N_{\text{bin }}+N_{\text{sing }}}$. The multiple fraction is defined as $$f_{\mathrm{m}\mathrm{u}\mathrm{l}}=\frac{N_{\mathrm{m}\mathrm{u}\mathrm{l}}}{N_{\mathrm{m}\mathrm{u}\mathrm{l}}+N_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g}}}.$$(10)

Here N_sing, N_bin, and N_mul are the number of single stars, binary systems, and multiple systems, respectively.

The initial distribution of the masses of the two binary components for lower-mass binaries is shown in Fig. 1. It is important to note that the binaries are first generated from a list of stellar masses drawn from the IMF and then assembled into a SC. Due to the dense environment inside the SC this can lead to an overlap of binary pairs or to multiple pairs of binaries being extremely close to each other. Therefore, not all of the stars originally set into the cluster as binaries turn out to be an actual binary pair as determined from their orbital energy. Consequently, the actual binary fraction is only around 90% in the initial timestep, f_bin = 0.9.

The clusters in all models are evolved over a period of ≈300 Myr by the code NBODY6 (Aarseth 1999, 2003), which couples regularising techniques (Kustaanheimo & Stiefel 1965; Aarseth & Zare 1974; Mikkola & Aarseth 1990) with a fourth-order Hermite predictor-corrector scheme (Makino 1991; Makino & Aarseth 1992). The code contains synthetic evolutionary tracks for both single star and binary star evolution (Tout et al. 1996; Hurley et al. 2000, 2002). The velocities for the SN kicks are drawn from a Gaussian distribution with σ = 40 km/s.

To distinguish between T1, which is caused by gas expulsion, and T2, caused by energy-equipartition-driven evaporation, Eq. (29) from Dinnbier & Kroupa (2020) was used to compute the threshold time, t_ee. The threshold time before which escaping stars are counted towards T1 is t_ee,main = 81.6 Myr for the main model. For the lower-mass models, we list the averages in Table 1. We note that t_ee slightly varies between individual models. For example for M2 t_ee ranges from 78.0 Myr to 86.8 Myr.

Fig. 1 Mass of the primary, m1, versus the mass of the secondary, m2. The grey lines are (from top to bottom) for the cases $q=\frac{m_2}{m_1}=1$, q = 0.6, and q = 0.3. The points show stars drawn randomly from the canonical IMF.

2.2 The computation of the observed masses

It is assumed that the observer is unable to resolve any of the binaries. Since the stars spend most of their time on the main sequence the MLR from Kroupa et al. (1993, in their Appendix A) can be used to get the luminosity, L(m), for each star (see also Reid et al. 2002). The main advantage of this MLR over others (e.g. Hertzsprung 1923; Russell et al. 1923; Eker et al. 2015) is that it extends down to very low-mass stars (0.07 M_⊙) and physically correctly accounts for the non-linearity near 0.35 M_⊙. The luminosity of the members of multiple systems is then added and converted back to a mass that represents the mass an observer would assume for the unresolved system: $$m_{\mathrm{o}\mathrm{b}\mathrm{s}}=L^{-1}(L\left(m_1\right)+L\left(m_2\right)),$$(11)

with the masses of the individual stars, m₁ and m₂, and the apparent mass an observer that cannot resolve the binary would compute from the total luminosity, m_obs.

Figure 2 shows the used MLR. The grey line follows the function M_V = (10.8 ± 0.1)(m/M_⊙)^{−0.491±0.006} + (−6.1 ± 0.1). This function was fitted using the data points for m > 0.96 M_⊙ and is used to compute the luminosities for all stars with m > 1.5 M_⊙. For all other masses the values of the closest red data point are used. The magnitudes, M_V, are transformed to luminosities, L, using the following equation: $$\frac{L(m)}{L_{\odot }}={10}^{-0.4(M_V(m)-M_{V,\odot })},$$(12)

where M_V,⊙ = 4.83 is the absolute V-band magnitude of the Sun.

It is important for our work to use an observationally derived MLR, rather than a theoretical one, since the observationally derived MLR correctly displays the KTG peak (Kroupa et al. 1990) of the stellar luminosity function at the transition from partially convective to fully convective stars (Kroupa 2002a). If not accounted for correctly, this leads to spurious results on the stellar MF (Kroupa et al. 1990).

In the following sections, the results for both the main and the lower-mass models are shown. The results for the lower-mass models are either averaged within their respective group (M2, M3, or M4; i.e. in the cases where quantities are shown over time) or the quantities of the lower-mass models are shown together (i.e. when the quantities of the individual binaries are shown).

Fig. 2 MLR used in this work. The absolute V-band magnitudes of stars, MV, are shown in relation to their stellar mass, m. The red dots show the data points from Kroupa et al. (1993), while the grey line shows the fit used for m > 1.5 M⊙.

3 Results and discussion

3.1 The apparent stellar masses of unresolved multiples

For the purpose of this discussion, all bright stars are assumed to be on the main sequence, that is to say, the short phases between when they are on the main sequence and when they become a black hole (BH), neutron star (NS), or white dwarf are ignored. Since the luminosity is added instead of the masses the relationship of the luminosity and the mass determines whether or not the mass of a system is over- or underestimated.

In Fig. 3, the L(m) used in this work and its first and second derivative are shown. We note that the anomaly in the second derivative at 1.5 M_⊙ stems from the transition from using the next data point from the MLR to using the fitted function. This leads to a small jump in the derivative of L(m), which is further enhanced through the second derivative. Since the second derivative is positive, L(m) must be a convex function in the given interval.

To understand, how this affects our estimate of the total mass of a binary, the following is considered: Assume two stars with masses m₁ and m₂. If the function is convex in an interval [0, m₁ + m₂] (i.e. every point lies below its chord), where L(0) = 0, it immediately follows that $$\frac{L(m_1+m_2)}{m_1+m_2}>\frac{L\left(m_1\right)}{m_1},$$(13)

Fig. 3 Luminosity of a star (Eq. (12) and Fig. 2) and its first and second derivative with respect to the stellar mass as a function of the stellar mass.

$$\frac{L(m_1+m_2)}{m_1+m_2}>\frac{L\left(m_2\right)}{m_2}.$$(14)

Multiplying with the denominator on the right side results in $$L(m_1+m_2)\frac{m_1}{m_1+m_2}>L\left(m_1\right),$$(15) $$L(m_1+m_2)\frac{m_2}{m_1+m_2}>L\left(m_2\right).$$(16)

Adding both equations together leads to $$L(m_1+m_2)\frac{m_1+m_2}{m_1+m_2}>L\left(m_1\right)+L\left(m_2\right),$$(17) $$L(m_1+m_2)>L\left(m_1\right)+L\left(m_2\right).$$(18)

Thus, the luminosities of the individual stars added up is always going to be smaller than the total luminosity of a single star with mass m₁ + m₂. This means that in this case the mass of an unresolved binary is always going to be underestimated. Analogous it could be shown that for a concave function with L(0) = 0 the mass of an unresolved binary is always going to be overestimated. An example of this can be seen in Fig. 4. It should be noted here that the L(m) used in the present work is only known down to 0.08 M_⊙. This can affect whether the mass of the binaries is over- or underestimated especially for the low-mass stars.

Figure 5 shows by how much an observer would underestimate the total mass of the multiple systems from the main model (top) and the models M2 (bottom) at the end of the simulation. The bottom plot shows all stars of the 5 models M2 together. The real mass of a system is m_rl = m₁ + m₂. The observed mass of the apparent star is lower for systems with higher q. The minimum of m_obs/m_rl for binaries with both components on the main sequence is ≈0.55. All systems below this minimum are either of an order higher than two or include BHs, NSs, or other non-main-sequence stars, which are assigned a luminosity of 0 L_⊙.

The underestimation of the masses of stellar systems leads to an underestimate of the total masses of the SC and its tails. In all three components (SC, T1, and T2) the fraction of higherorder multiples (more than two stars) is small. In fact, in the main model, no higher-order multiple composed of only main-sequence stars is found in the last timestep. While the number of binaries in the SC and its tails differ, the qualitative distribution of binary masses is similar in all three components.

Fig. 4 Examples on how the non-linearity of L(m) can lead to an underestimate or an overestimate of the total mass of a binary star system.

Fig. 5 Masses of multiple star systems in the SC (left), T1 (centre), and T2 (right) of the main model (top) and the five M2 models (bottom) as seen by an observer over the real total mass of the systems at the end of the simulation (T = 300 Myr). Binaries composed of only main-sequence stars are shown as dots, while multiple systems containing non-main-sequence stars are displayed as crosses. The colour denotes the mass ratio of the system.

3.2 Evolution of the binary fraction and the total photometric and dynamical mass

Figure 6 shows the number of stars in single and multiple systems for the SC, T1, and T2 (left panel) and the fraction, f_mul, of multiple systems over the total number of systems (right panel). A system can be a single star or a multiple system. We note that a vast majority of the multiple systems are binaries and only a few of them are triple and quadruple systems (see also Fig. 5). For correctness, they will still be referred to as multiple systems in this work. The number of stars in multiple systems (initially close to 100%; see also Fig. A.1 for a linear plot that contains T = 0) decreases rapidly during the first megayear of cluster evolution, which is due to dynamical encounters in the dense environment (see also Kroupa 1995a,b). There are a few dynamical ejections from the SC before gas expulsion starts at 0.6 Myr.

With gas expulsion T1 forms rapidly and inherits the fraction of stars in multiples from the young SC models. The low multiple fraction in T1 at the beginning of its formation is due to dynamical ejections as stars ejected dynamically before t_ee are counted towards T1. In the main model the higher density in the SC leads to a decrease in the multiple fraction, while it remains constant at f_mul ≈ 0.38 in T1. T2 evaporates from the SC more slowly and, therefore, at the end of the simulation, has almost the same fraction of multiples (f_mul ≈ 0.31) as the ageing SC (f_mul ≈ 0.32). During the formation of T2 the multiple fraction fluctuates due to low-number statistics. Since the simulated SC moves in the Galactic plane, the tails form in the Galactic plane as well.

Figure A.2 shows the models M3 (top) and M4 (bottom). Due to their lower densities and longer crossing times the binaries are dissolved more slowly and to a lesser extent the less massive a SC becomes (Kroupa 1995a,b). This leads to higher multiple fractions in T1 and the SC. It is also noticeable that the difference between T1 and the SC is decreased. This is because the time after T1 formed is not sufficient to reduce the multiple fraction in the SC further. It is also worth noting that in the models M3 and M4 the binary fraction in the SC increases towards the end of the simulation. This is due to the preferential evaporation of single stars. The overall binary fraction in these models is decreasing up to the end.

However, the lower-mass models lose a larger fraction of their mass to T1, so that at the end of the simulation the number of T1 stars is larger than the number of SC stars. This is consistent with their shorter relaxation times (see Table 2). This is caused by the less deep gravitational potential counteracting the tidal forces on the SC.

Figure 7 shows the development of the masses of the SC, M_SC, and the tidal tails, M_T1 and M_T2, of the main model. The total mass (shown as the thick solid line) includes all stars, NSs, and BHs. For the stars and observed masses only stars on the main sequence are taken into account. This simplification is possible due to the short giant phases and low mass-fraction of high-mass stars in the SC and its tails. The mass white dwarfs would add leads to no noticeable change. Therefore, white dwarfs are neglected as well. The dashed curve takes into account that an observer would not only underestimate the SC mass due to objects they cannot see, but also due to unresolved multiples. As explained in Sect. 3.1, the apparent mass of a multiple system is smaller than the sum of the masses of the companions, resulting in a lower ‘measured’ total mass. Therefore, this curve only includes stars on the main sequence and in the case of multiples the apparent mass of the star an observer would deduce from the luminosity of the unresolved system. The exact values of the underestimate at the end of the simulations can be found in Table 2.

The sum of the total masses of the SC and its tails $(M_{\mathrm{S}\mathrm{C}}^{\text{total }}+M_{\mathrm{T}1}^{\text{total }}+M_{\mathrm{T}2}^{\text{total }})$ decreases over time due to mass loss through stellar evolution. The SC loses mass continuously through ejection, evaporation and, to a lesser extent, stellar evolution. T1, on the other hand, initially gains a lot of mass and then slowly loses mass due to stellar evolution only. This is because gas expulsion only plays a role in the early stages of the evolution (see also Fig. 6), after which mass is lost due to stellar evolution. T2 gains mass almost continuously, which means that evaporation of stars from the SC adds mass faster than can be lost due to stellar evolution.

In addition to the underestimate of the SC, T1, and T2 masses due to the above-described bias from unresolved multiple systems, the amount by which the mass of a SC would be overestimated using the velocity dispersion is computed. As mentioned in the introduction, Illingworth (1976) and Spitzer (1987) found that the determined dynamical mass is proportional to the square of the velocity dispersion of the centres of mass of the stellar systems in the SC, σ_COM. Following Rastello et al. (2020) it is assumed that an observer would see the centre of light movement, rather than the centre of mass movement. The proportionality factor between the cluster mass and the velocity dispersion depends on the assumed mass profile of the cluster. To identify merely the error caused by the velocity dispersion, we therefore used $$M_{\mathrm{S}\mathrm{C}}^{\text{dynamical }}=M_{\mathrm{S}\mathrm{C}}^{\text{total }}{\left(\frac{{\sigma }_{\mathrm{C}\mathrm{O}\mathrm{L}}}{{\sigma }_{\mathrm{C}\mathrm{O}\mathrm{M}}}\right)}^2,$$(19)

with σ_COL being the velocity dispersion of the centre of lights of the stellar systems. For our calculations we use the velocity dispersions in z-direction (σ_COL = σ_COL,z, σ_COL,z, σ_COM = σ_COM,z).

Here it is important to determine which stars an observer would count as members of the SC, as stars with too large radial velocities are often discarded as non-members. To this end, the stars are sorted by their distance to the SC centre. For each star the velocity dispersion σ_COL,z is computed from its own velocity and the velocities of its 20 neighbours in each direction. The velocity dispersion is computed for each star individually to account for the velocity dispersion being larger in the centre of the SC, than in the outskirts. Systems whose velocities are larger than 3σ_COL,z are excluded from the calculations. This method is applied iteratively until no stars remain outside this variance.

As seen in Fig. 7, the dynamical mass is significantly larger than the total mass of the SC. For the main model, we measure an overestimate by a factor of 1.86 (Table 2) at the end of the simulation.

At the same time, the mass of the SC is underestimated by 3 % due to the presence of unseen BHs and NSs and a further 15 % due to unresolved binaries. The effect caused by BHs and NSs is small since many of them get ejected from the SC due to SN kicks. Similarly, the mass in T1 (T2) is underestimated by 14 % (6 %) due to not detecting BHs and NSs and a further 14 % (14 %) due to unresolved binaries. Here, the effect of BHs and NSs is stronger, since the remnants kicked out from the SC by SNe are added to its tails in our model. However, since these remnants are kicked out with high velocities (≈40 km/s), their phase space positions do not necessarily match those of the other tail members. Since the most massive stars are short-lived and the majority of SNe explode early on, it is expected that dark objects play a more significant role in T1’s mass underestimate, while in the SC and T2 unresolved binaries cause the majority of the underestimate.

Figure A.3 shows a similar development for the Model M4. Due to the higher binary fraction the dynamical overestimate of the SC mass is higher than for the main model (17.52 compared to 1.86 for the main model; Table 2). The underestimate of the SC caused by BHs, NSs, and binary systems is similar to that of the main model at 22%. We also note that the dynamical mass of the SC increases at the end of the simulation together with the binary fraction.

Figure 8 shows the under- and overestimates of the masses at different times for all models used in this work. These values found for the overestimates using the dynamical mass estimate are smaller than the results found by Rastello et al. (2020), who for a cluster of a similar binary fraction as the ones shown in this work (M_SC,tot = 275 M_⊙, f_mul = 0.421) found that the velocity dispersion is overestimated by a factor of 6.7. Due to the square in Eq. (19), this would correspond to a mass overestimate by a factor of 44.8. The main reason for that difference is that in the present work the stars with the largest velocities are removed as explained above, since this work assumes that an observer would probably not count these stars as cluster stars. This leads to a smaller velocity dispersion in the present work. A full list of the over- and underestimates of the SC masses at the end of the simulations can be found in Table 2.

Both the overestimate of the dynamical mass and the underestimate of the visible mass of the SCs become more pronounced over time. It is visible that the dynamical mass is both mass and time dependent. This is mainly due to the high initial cluster density leading to a large velocity dispersion. If a SC already has a large velocity dispersion, binaries are not increasing the perceived velocity dispersion by a large factor. Low-density, low-mass clusters on the other hand have a low velocity dispersion, which means that binaries increase the perceived velocity dispersion and the dynamical mass by a large factor.

As mentioned in the introduction, observers often use the line-of-sight velocity dispersion to determine the dynamical mass instead of the photometric mass. This work clearly demonstrates that not only is the dynamical mass overestimated due to binaries; the photometric mass is significantly underestimated. This increases the apparent discrepancy between the cluster masses determined by different methods. This is relevant for the Hyades star cluster, which is super-virial by a factor of about two (Kroupa et al. 2022, and references therein). For the main model of this work, the dynamical and luminous mass differ initially by a factor of 1.6, which increases to 2.3 at the end of the simulation. For the other models this factor is even larger. However, many of the binaries in Hyades have been resolved (Gizis & Reid 1995), leaving a smaller difference than in our models that assume all binaries are unresolved.

Fig. 6 Numbers of single stars and stars in multiple systems for the SC, T1, and T2 (left). The fraction of multiple systems over the total number of systems (counting single stars as systems), fmul (Eq. (10)), computed from these numbers can be seen on the right. The top panels show the main model, while the bottom panels are for the five M2 models (Table 1).

Fig. 7 Development of the total mass, mass in main-sequence stars, and observed mass of the SC (top), T1 (middle), and T2 (bottom) of the main model over time. The dynamical mass estimate according to Eq. (19) is only shown in the top panel for the SC.

Fig. 8 Over- and underestimate of the SC masses due to binary fractions and dark objects at different times depending on the cluster’s current mass. The values are averaged over the last five timesteps, which span less than 1.3 Myr in total.

3.3 The velocity dispersion of Galactic field stars

By analysing star-count data, Li et al. (2023) deduced apparently systematic variations in the IMF with the velocity dispersion in z-direction (perpendicular to the Galactic plane) and metallicity. Such a variation in the IMF has important astrophysical implications (Kroupa et al. 2013). Here we test if the origin of stars in star clusters may not partially or wholly account for this apparent variation. This may result if the binary fraction and velocity dispersion conspire to lead to apparently systematic variations in the observationally deduced IMF.

The velocity dispersion vertical to the Galactic plane is used as a proxy for the line-of-sight velocity dispersion. This is done because the velocity ellipsoid in radial direction with respect to the Galactic centre is larger than in tangential direction. However, the orientation of these axis in the non-rotating x-y plane changes over time. Therefore, the velocity dispersion in x- and y-direction is expected to vary with time, as the population revolves around the Galaxy (for a more detailed description of how the tails develop over time, see Dinnbier & Kroupa 2020). In the case of multiples the velocity of the centre of mass is used for the calculations. We note that the observer would observe the SC primarily from within the Galactic disc and, therefore, would have to account for the positions of the SC relative to the Sun.

Snapshots of the first eight timesteps of the main model are shown in Fig. 9. The tidal radius, calculated according to Eq. (1) from Baumgardt & Makino (2003) is shown as a black circle. Only one star escapes the SC at 0.3 Myr, which is clearly dynamically ejected. However, in the same timestep we also see a star within the cluster with a velocity exceeding 100 km/s. This is typical for an ejected star about to leave. After 0.6 Myr it is difficult to distinguish between stars that left the SC primarily due to gas expulsion and those leaving primarily due to dynamical effects. It is clear, however, that only a handful of stars leave the SC before gas expulsion. These stars leaving the SC propagate with a velocity of a few tens to hundreds of km/s.

The velocity dispersion for stellar systems in the SC and its tails is shown in Fig. 10. The high peak in the velocity dispersion at the very beginning for both the SC and T1 is caused by stars being ejected due to stellar-dynamical encounters. As gas expulsion starts the stars lost due to gas expulsion quickly form the majority of the T1 stars. These stars oscillate around the Galactic plane. Since the T1 stars are all expelled during a very short time interval through gas expulsion, this translates to oscillations of the velocity dispersion for T1 around about 3 km/s. They can be distinguished easily from the stars in the SC and T2 due to their larger velocity dispersion.

The T1 stars complete about 9.9 vertical oscillations over the 300 Myr of simulation. This translates to a period of about 32 Myr, which is comparable to the vertical period of the Sun around the Galactic plane (Bahcall & Bahcall 1985). T1 ends up far more dispersed in both the Galactic plane and the orthogonal direction, as can be seen from the velocity dispersions. Therefore, T2 would be expected to be much easier to observe than T1.

As can be seen in Table 2 differences between the velocity dispersions of different models are small compared to the differences between T1 and T2. More massive SCs produce tidal tails with fewer binaries and a larger velocity dispersion. Only looking at this result, it might seem as if field populations with a higher velocity dispersion would automatically be expected to host fewer binaries. However, we have also seen that T1 hosts a larger fraction of binaries and has a larger velocity dispersion than T2. For the models investigated here, the difference between different tails is larger than the difference produced by different SC masses. To get a complete quantification for a correlation of binary fraction and velocity dispersion in the Galactic field a complete simulation of SCs following the MF of embedded star clusters is needed. However, it should be noted that Li et al. (2023) do find a higher binary fraction in dynamically hot stellar populations than in cold populations observed using LAMOST. They distinguish the two groups of stars by their action in the z-direction, J_z, and call those with J_z > 20 kpc km/s dynamically hot. Stars with line-of-sight velocities exceeding 25 km/s are assumed to be in binaries by them.

Combining the previous results – that the binary fraction in T1 is larger than in T2 and the SC – with the fact that the stars in T1 have the largest σ_z, it follows that the binary fraction in groups of stars belonging to the same tail is expected to be larger if σ_z is large, and vice versa. In Fig. 11 the fraction of multiple systems, f_mul, is shown depending on the velocity dispersion. To obtain this plot several cylindrical volumes with a radius of 250 pc and infinite length, oriented orthogonal to the Galactic plane, are inserted into the tidal tails at different angular distance from the SC as seen from the Galactic centre. These cylinders overlap. The reason for using cylinders is that the velocity dispersion in x- and y-direction (in the non-rotating frame) is affected by the Galactic rotation, while in z-direction it is independent of it. Within each of those cylinders the fraction of multiple stars, f_mul, and σ_z is measured. For each timestep, 100 cylinders are used; however, they are ignored if they contained fewer than 100 systems. From the colour coding it is visible that both σ_z and f_mul are increasing with increasing distance to the SC. However, the spread also increases significantly, so that it would become difficult for an observer to detect a correlation between σ_z and f_mul at large distances to the SC. Especially in the trailing tail, the number of binaries seems to decrease again at large distances (10°) to the SC. However, this decrease corresponds to about 10 pairs per cylinder. Because of their size, the cylinders overlap, so that the decrease in the binary fraction is likely caused by a local under-density of binaries in this particular model snapshot.

The velocity dispersion in this snapshot (shown in Fig. 11) is larger in the trailing tail, than it is in the leading tail. However, as already shown in Fig. 10 the velocity dispersion in T1 fluctuates over time. In Fig. 12, the same plot as Fig. 11 is shown for T = 283 Myr. Here the leading tail has a larger velocity dispersion than the trailing one. The phase difference is caused by the leading tail initially moving towards the Galactic centre and the trailing tail moving away from it.

Overall these results show that there is a correlation between binary fraction and σ_z. However, a stellar-dynamical population model of the Galactic field population, consisting of many accumulated populations, would need to be computed to estimate systematic biases in deriving the stellar IMF from star counts.

Fig. 9 Snapshots of the SC (main model; Table 1) in the x-y plane (the Galactic centre being in negative x-direction) for the first eight timesteps to show the expansion due to gas expulsion. The tidal radius is shown as a black circle, and the stars are colour-coded according to the absolute value of their velocities relative to the SC’s centre of mass.

Fig. 10 Velocity dispersion of stars in the z-direction for all stars between 0.3 and 0.7 M⊙ in the SC (blue), T1 (orange), and T2 (green) for the main model.

Fig. 11 Fraction of stars in multiple systems over the velocity dispersion of stars in the z-direction at the end of the simulation (T = 300 Myr) for the main model. The points are colour-coded according to the angular distance in degrees of the cylinder from the SC with respect to the Galactic centre. Positive values of ϕSC indicate the stellar populations in the leading T1 and T2 tails, while ϕSC < 0 indicate the trailing tail ones.

Fig. 12 Same as Fig. 11, but at T = 283 Myr.

3.4 The effect of the binary fraction on the stellar MF

The stellar MF is fitted to the stars in the generic mass range between 0.3 and 0.7 M_⊙. This is the typical mass range used by other researchers (see Table 2 in Yan et al. 2024; Paresce & De Marchi 2000; Baumgardt & Hilker 2018; Lardo et al. 2022). The low-mass end is limited by completeness, while the highmass end serves to remove stellar evolution issues. For the entire simulated period of 300 Myr, these stars never leave the main sequence. Therefore, no post-main-sequence objects need to be taken into account and the MLR from Kroupa et al. (1993) can be used. In this work, two main methods for deriving the power-law index, α, were employed:

The first method is histogram fitting, which was mentioned by Maschberger & Kroupa (2009). It is visualised in Fig. 13. The stars are first sorted into logarithmic mass bins and their number is divided by the bin width, w = log₁₀(m_j+0.5) − log₁₀(m_j−0.5) (m is in units of M_⊙). j is the number of the bin and w is equal for all bins. The mass MF dN = ξ(m)dm = km^−αdm is a power-law, which becomes linear in log-log. Therefore, we translated the MF into the double logarithmic frame and obtain for each bar in the histogram at the position log₁₀(m) of height N_j: $$\xi ({\log }_{10}(m))=k{10}^{(1-\alpha ){\log }_{10}(m)}\ln (10),$$(20) $$N_j={\displaystyle \underset{{\log }_{10}\left(m_{j-0.5}\right)}{\overset{{\log }_{10}\left(m_{j+0.5}\right)}{\int }}}\xi ({\log }_{10}(m))d{\log }_{10}(m)$$(21) $$=\frac{k}{1-\alpha }({10}^{(1-\alpha )\log \left(m_{j+0.5}\right)}-{10}^{(1-\alpha )\log \left(m_{j-0.5}\right)}),$$(22) $$=\frac{k}{1-\alpha }{10}^{(1-\alpha ){\log }_{10}\left(m_j\right)}[{10}^{(1-\alpha )w/2}-{10}^{-(1-\alpha )w/2}],$$(23) $${\log }_{10}(N_j)=(1-\alpha ){\log }_{10}(m_j)+const.$$(24)

The function in Eq. (24) is then fitted to the data to obtain the MF.

The second method is the maximum likelihood estimator (MLE) described in Appendix A of Khalaj & Baumgardt (2013). The MLE is based on the work of Clauset et al. (2009); Maschberger & Kroupa (2009) who developed it for cases with only a lower boundary. Their main idea is to formulate the likelihood that a certain distribution of data points results from a MF with a given α and then maximising the function for α. As a first step both the fitting method and the MLE are tested on an optimally sampled canonical IMF (Kroupa et al. 2013; Yan et al. 2017) between m = 0.3 and 0.7 M_⊙. This gives α = 1.65 ± 0.03 for the fitting approach above and α = 1.62 ± 0.03 using the MLE. It is important to note here that there is no ‘correct’ value for the power-law slope here, since the true slope is a two-part power-law instead of the one-part power-law that is fitted here. The lower value from the MLE stems from the MLE putting more emphasis on matching the low-mass portion of the MF than fitting a histogram does. If the low-mass portion of the MF would have a larger α than the high-mass one we find that α determined with the MLE to be larger than the one deduced from fitting a histogram.

Figure 14 shows how the power-law index of the MF for stars between 0.3 and 0.7 M_⊙, ${\alpha }_{0.3-0.7M_{\odot }}$, develops over time. In all 3 groups (SC, T1 and T2) the observed MF is flatter than the real one that would be measured if all binaries were resolved by about Δα = 0.2. For the main model the real MF for stars in the SC is steeper than the one for T1 (${\alpha }_{0.3-0.7M_{\odot }}$ = 1.8 v.s.1.5). For the SC and T1 the real MF stays constant for the majority of the simulation. Only in the very beginning, when T1 starts forming, are variations visible. However, the MF of the SC is not always steeper than that of T1. Looking at the lower-mass models, we find no systematic tendency in either direction. The observed variations can be attributed to the MF not being optimally, but rather stochastically, sampled in the models used here and the models being non-mass-segregated (Kroupa et al. 2013). Compared to the MF of the SC and T1, the MF of T2 varies significantly. This is due to the slow evaporation process that forms T2. As more stars join T2, the variations decrease and the MF of T2 becomes comparable to the other two.

The lower panels of Fig. 14 show the α-values for the same data as the upper panels, but obtained using the MLE. In all cases the real and observed values for α are closer together. It is also visible that α of the SC and T1 remains below the value obtained for an optimally sampled IMF, despite this not being the case if histograms are used to determine α. This is a result of the clusters being stochastically rather than optimally sampled. These differences in the direction in which the α value of a stochastically sampled MF deviates from the optimally sampled one is the result of the MLE putting more emphasis on the low-mass portion of the MF.

Figure 15 shows how the fraction of stars in multiples affects the difference between the real α of the population and the α an observer who cannot resolve the binaries measures, Δα = α_rl − α_obs for the cylinders introduced in Sect. 3.3. The values of α used here were obtained via fitting a histogram (top) or using the MLE (bottom). Using a weighted least square fit, a linear function, $\mathrm{\Delta }{\alpha }_{0.3-0.7M_{\odot }}=b{f}_{\mathrm{m}\mathrm{u}\mathrm{l}}$, was fitted.χ² for a linear function, y = f(x) = bx is $${\chi }^2={\displaystyle \sum _{i=1}^N}\frac{{(y_i-b{x}_i)}^2}{{\sigma }_i^2},$$(25)

with the fitting parameter, b, and the error for y_i, σ_i. From this follows $$b=\frac{{\displaystyle \sum _{i=1}^N}\frac{x_i{y}_i}{{\sigma }_i^2}}{{\displaystyle \sum _{i=1}^N}\frac{x_i^2}{{\sigma }_i^2}}\text{ and }\mathrm{\Delta }b=\frac{1}{\sqrt{{\displaystyle \sum _{i=1}^N}{\left(\frac{x_i}{{\sigma }_i}\right)}^2}}.$$(26)

Using f_mul as x and $\mathrm{\Delta }{\alpha }_{0.3-0.7M_{\odot }}$ as y, this yields b = 0.8 ± 0.3 for the main model. The values for all models are listed in Table 2. The positive correlation matches the results by Li et al. (2023) qualitatively, but not quantitatively. This is due to largely different methods being used in their work, compared to the present one. In addition to using a mass range between 0.1 and 3.0 M_⊙, compared to the 0.3 to 0.7 M_⊙ used in the present work, Li et al. (2023) neglected the luminosity contributed by the secondary instead of adding the luminosities.

The colour coding in the left panel of Fig. 15 shows that the multiplicity and $\mathrm{\Delta }{\alpha }_{0.3-0.7M_{\odot }}$ are larger for cylinders with a larger real α. This seems counter-intuitive, since T1 has a flatter observed MF and higher multiplicity than T2 (see Figs. 6 and 14). However, looking at the right panel it becomes clear that the internal spread within T1 dominates over the differences between T1 and T2. Due to the stars in T1 being spread out more than the stars in T2, the number of stars per cylinder is much lower. This is reflected in the large error bars of the cylinders containing mainly T1 stars.

Fig. 13 Bins used to fit the MF between 0.3 and 0.7 M⊙ to the stars of the SC in the initial timestep (T = 0) of the main model. The fitted MF is shown with a solid line, and the MF determined using the MLE (with k fitted to best fit the bins) with a dashed line. The error bars show the Poisson uncertainties.

Fig. 14 MF slope obtained using the fitting approach (Eq. (24), top) and the MLE (bottom) between 0.3 and 0.7 M⊙ over time for the main model. Both the MFs as seen by an observer who cannot resolve the binaries (dashed) and the actual real MFs (solid) of the resolved SC, T1, and T2 are shown. The coloured areas show the 1σ confidence span. The canonical value of 1.65 is shown as a dotted black line.

Fig. 15 Difference between the observed and real α between 0.3 and 0.7 M⊙ depending on the multiplicity for the last timestep of the main model (top) and the low-mass models (bottom). The data is colour-coded according to the real α (left) and the portion of T1 stars (right). Populations with a larger fraction of unresolved multiple systems thus have a larger Δα and therefore a smaller observed α, i.e. the observer deduces a flatter MF. The grey area shows the uncertainty in the slope. Note that for the five low-mass models, all data points are plotted and fitted together.

4 Conclusions

In the present work, realistic N-body simulations of SCs with masses between 800 M_⊙ and 6400 M_⊙ carried out with NBODY6 were used to determine the effect of multiple star systems on the measured masses and MFs. These masses are near the upper limit of masses of nearby SCs (Ebrahimi et al. 2022). The apparent mass a multiple system would have was computed by adding the luminosities produced by the companion stars and using the MLR from Kroupa et al. (1993, see also Reid et al. 2002). The main results are as follows:

(1) Stars in T1 have a higher fraction of multiples than those in T2 and the SC. This is due to the T1 stars being expelled early, before the multiple fraction in the SC drops to its final value.

(2) Due to the presence of remnants (NSs and BHs) and unresolved multiples, the photometric mass underestimates the total mass of the SC by nearly one-fifth. Binaries are the main contributor to this effect. In the main model, they cause a 14–15% underestimate of the total photometric mass of the SC and its tails. BHs and NSs have only a small effect in the SC, where they cause a 3% underestimate, but a very strong effect in T1, the total photometric mass of which is underestimated by an additional 14% due to dark objects. The dynamical mass of the SC in the main model is a factor of 1.6 higher than the luminous mass at the beginning of the simulation, and a factor of 2.3 higher at its end. This is reminiscent of the Hyadesmass problem discussed in Röser et al. (2011) and Kroupa et al. (2022).

(3) The MF in the stellar mass range 0.3–0.7 M_⊙ appears more top-heavy to an observer who cannot resolve multiple systems. In both tails and the SC, the difference $\mathrm{\Delta }{\alpha }_{0.3-0.7M_{\odot }}$ is 0.2 ± 0.3. While this is within 1σ, the fact that this difference is present at all timesteps, in both the SC and its tails and in all simulations, suggests a systematic effect.

(4) While there is a positive correlation between f_mul, σ_z, and $\mathrm{\Delta }{\alpha }_{0.3-0.7M_{\odot }}$ within a single model, when looking at different models one can find a positive correlation between the initial SC mass and the velocity dispersion in the tidal tails. However, the correlation between the SC mass and binary fraction is negative.

In this work it is assumed that all binaries are unresolved; however, in reality some binaries might be resolved and some not. Therefore, future studies should take the distance of the SC to the observer into account when studying these effects.

The field stars we observe in the Galaxy consist of a mixture of tidal tails of SCs of many different initial masses. To recreate the stellar population actually observed in the field, a full population synthesis that includes a range of models with different SC masses would be necessary (see e.g. Kroupa 2002b; Marks & Kroupa 2011; Dabringhausen et al. 2022). This would have to be done by including a range of IMFs that reflect their dependence on the metallicity and density of the gas cloud the SC is born from, as was found by comparing the results of N-body models to observed SC properties, such as their masses, half-mass radii, and binary fractions (Marks et al. 2012), and by comparing the chemical evolution of galaxies computed for different IMFs to the metallicities observed in dwarf galaxies (Lee et al. 2009; Gunawardhana et al. 2011; Yan et al. 2020). The two methods lead to similar results for the upper end of the IMF. It is also important to mention that Baumgardt et al. (2023) measured the upper end of the IMF of 14 Small and Large Magellanic Cloud clusters directly. While the masses and metallicities of these clusters should theoretically lead to a flattened IMF, their errors are too large to detect any systematic variation. Similarly, Dickson et al. (2023) investigated 37 Galactic globular clusters without finding any evidence for a variation in the IMF after fitting multi-mass models to them. However, they only cite statistical errors and, in contrast to the analysis by Marks et al. (2012), Dickson et al. (2023) neglect very early stellar dynamical processes during and before gas expulsion, which makes it difficult to assess the actual errors of the found MF slopes. Gas expulsion itself is not taken into account in these models.

The work presented here thus points towards the required deep analysis needing to be performed when studying the field-star population in order to infer how the stellar IMF varies with the physical conditions of star formation.

Acknowledgements

We thank an anonymous referee for many helpful suggestions. The authors acknowledge support from the Grant Agency of the Czech Republic under grant number 20-21855S and through the DAAD-Eastern-Europe Exchange grant at Bonn University. HW thanks Roberto Capuzzo Dolcetta for valuable discussions.

Appendix A Supplementary figures

Fig. A.1 Same as Fig. 6, but with a linear timescale.

Fig. A.2 Same as Fig. 6, but for the models M3 (top) and M4 (bottom).

Fig. A.3 Same as Fig. 7, but the average for the 16 M4 models (Table 1).

References

All Tables

All Figures

	Fig. 1 Mass of the primary, m1, versus the mass of the secondary, m2. The grey lines are (from top to bottom) for the cases $q=\frac{m_2}{m_1}=1$, q = 0.6, and q = 0.3. The points show stars drawn randomly from the canonical IMF.
In the text

	Fig. 2 MLR used in this work. The absolute V-band magnitudes of stars, MV, are shown in relation to their stellar mass, m. The red dots show the data points from Kroupa et al. (1993), while the grey line shows the fit used for m > 1.5 M⊙.
In the text

	Fig. 3 Luminosity of a star (Eq. (12) and Fig. 2) and its first and second derivative with respect to the stellar mass as a function of the stellar mass.
In the text

	Fig. 4 Examples on how the non-linearity of L(m) can lead to an underestimate or an overestimate of the total mass of a binary star system.
In the text

Fig. 5 Masses of multiple star systems in the SC (left), T1 (centre), and T2 (right) of the main model (top) and the five M2 models (bottom) as seen by an observer over the real total mass of the systems at the end of the simulation (T = 300 Myr). Binaries composed of only main-sequence stars are shown as dots, while multiple systems containing non-main-sequence stars are displayed as crosses. The colour denotes the mass ratio of the system.

In the text

	Fig. 6 Numbers of single stars and stars in multiple systems for the SC, T1, and T2 (left). The fraction of multiple systems over the total number of systems (counting single stars as systems), fmul (Eq. (10)), computed from these numbers can be seen on the right. The top panels show the main model, while the bottom panels are for the five M2 models (Table 1).
In the text

	Fig. 7 Development of the total mass, mass in main-sequence stars, and observed mass of the SC (top), T1 (middle), and T2 (bottom) of the main model over time. The dynamical mass estimate according to Eq. (19) is only shown in the top panel for the SC.
In the text

	Fig. 8 Over- and underestimate of the SC masses due to binary fractions and dark objects at different times depending on the cluster’s current mass. The values are averaged over the last five timesteps, which span less than 1.3 Myr in total.
In the text

	Fig. 9 Snapshots of the SC (main model; Table 1) in the x-y plane (the Galactic centre being in negative x-direction) for the first eight timesteps to show the expansion due to gas expulsion. The tidal radius is shown as a black circle, and the stars are colour-coded according to the absolute value of their velocities relative to the SC’s centre of mass.
In the text

	Fig. 10 Velocity dispersion of stars in the z-direction for all stars between 0.3 and 0.7 M⊙ in the SC (blue), T1 (orange), and T2 (green) for the main model.
In the text

Fig. 11 Fraction of stars in multiple systems over the velocity dispersion of stars in the z-direction at the end of the simulation (T = 300 Myr) for the main model. The points are colour-coded according to the angular distance in degrees of the cylinder from the SC with respect to the Galactic centre. Positive values of ϕSC indicate the stellar populations in the leading T1 and T2 tails, while ϕSC < 0 indicate the trailing tail ones.

In the text

	Fig. 12 Same as Fig. 11, but at T = 283 Myr.
In the text

	Fig. 13 Bins used to fit the MF between 0.3 and 0.7 M⊙ to the stars of the SC in the initial timestep (T = 0) of the main model. The fitted MF is shown with a solid line, and the MF determined using the MLE (with k fitted to best fit the bins) with a dashed line. The error bars show the Poisson uncertainties.
In the text

Fig. 14 MF slope obtained using the fitting approach (Eq. (24), top) and the MLE (bottom) between 0.3 and 0.7 M⊙ over time for the main model. Both the MFs as seen by an observer who cannot resolve the binaries (dashed) and the actual real MFs (solid) of the resolved SC, T1, and T2 are shown. The coloured areas show the 1σ confidence span. The canonical value of 1.65 is shown as a dotted black line.

In the text

Fig. 15 Difference between the observed and real α between 0.3 and 0.7 M⊙ depending on the multiplicity for the last timestep of the main model (top) and the low-mass models (bottom). The data is colour-coded according to the real α (left) and the portion of T1 stars (right). Populations with a larger fraction of unresolved multiple systems thus have a larger Δα and therefore a smaller observed α, i.e. the observer deduces a flatter MF. The grey area shows the uncertainty in the slope. Note that for the five low-mass models, all data points are plotted and fitted together.

In the text

	Fig. A.1 Same as Fig. 6, but with a linear timescale.
In the text

	Fig. A.2 Same as Fig. 6, but for the models M3 (top) and M4 (bottom).
In the text

	Fig. A.3 Same as Fig. 7, but the average for the 16 M4 models (Table 1).
In the text