© The Authors 2025

1 Introduction

Globular clusters (GCs) were once considered the best approximations of simple stellar populations (e.g., Renzini & Buzzoni 1986). However, various studies have revealed groups of GC stars with different chemical compositions. It is now widely accepted that most GCs exhibit multiple stellar sequences and complex horizontal-branch morphology in color-magnitude diagrams (CMDs) built with appropriate photometric bands. All these features are nothing but different facets of the same phenomenon, known as the multiple stellar population phenomenon, according to which GCs host two main stellar populations: a first population (1P) with a Galactic-field-like chemical composition and a second population (2P) enhanced in He, N, Na, and depleted in C and O (see reviews by Kraft 1994; Gratton et al. 2012; Bastian & Lardo 2018; Milone & Marino 2022).

The formation of multiple stellar populations is a topic of intense debate. According to most scenarios, 2P stars formed in the central regions of the primordial GCs, whereas 1P stars had a sparser initial spatial distribution (see, e.g., D’Ercole et al. 2008; Bekki 2010, 2011; Renzini et al. 2015; Calura et al. 2019; McKenzie & Bekki 2021; Lacchin et al. 2022). In this context, binary stars can be valuable tools for constraining the spatial properties of 1P and 2P stars at the time of their formation and during their subsequent dynamical evolution. For instance, in high-density environments, frequent stellar encounters lead to a high disruption rate of binary systems. Conversely, in low-density environments, the binary fraction is expected to be higher. Therefore, the incidence of binaries among 1P and 2P stars can provide insight into the properties of the environments where these stellar populations formed (e.g., Vesperini et al. 2011; Hong et al. 2015, 2016; Sollima et al. 2022; Hypki et al. 2022).

Early determinations of the incidence of binaries among multiple populations were based on spectroscopy. Lucatello et al. (2015) analyzed 21 binaries located in 10 GCs around the half-light radius (R_hl) and determined that the fraction of 1P binaries is 4.1 ± 1.7 times higher than that of 2P binaries. Based on 12 binaries in NGC 6362 located beyond 0.5 R_hl, Dalessandro et al. (2018) found their fraction among 1P and 2P stars to be 4.7 ± 1.4 and 0.7 ± 0.7%, respectively. Similar conclusions were obtained by D’Orazi et al. (2010), who found that four out of the five studied barium stars of NGC 6121 are sodium-poor and oxygen-rich.

Conversely, Milone et al. (2020) investigated the innermost regions of five GCs, namely NGC 288, NGC 6121, NGC 6352, NGC 6362, and NGC 6838 from 0.3 R_hl to about 0.8 R_hl. Except for NGC 6121, the inner regions were found to be characterized by similar 1P and 2P binary fractions. Observations also provided the first evidence of the existence of mixed binaries composed of one 1P star and one 2P star, a population of binaries predicted by Hong et al. (2015, 2016) to form during exchange encounters in binary-single and binary-binary interactions.

In this paper we analyze the binary systems among the multiple populations of NGC 288, a GC that has been extensively studied in the context of multiple stellar populations. Previous studies based on both spectroscopy and photometry have identified two groups of stars with different abundances of certain light elements, including helium, carbon, nitrogen, oxygen, sodium, and aluminum (e.g., Shetrone & Keane 2000; Kayser et al. 2008; Smith & Langland-Shula 2009; Carretta et al. 2009; Milone et al. 2018; Marino et al. 2019). Additionally, photometric studies have revealed two distinct sequences in the CMD that can be traced continuously across various evolutionary stages, including the red-giant branch, subgiant branch, and main sequence (MS) (Lee et al. 2009; Piotto et al. 2013; Monelli et al. 2013; Milone et al. 2017a; Lagioia et al. 2021; Jang et al. 2022; Dondoglio et al. 2022).

All previous studies on binaries analyzed multiple systems within 1 R_hl and with a primary star mass $(M_{bin}^1)$ higher than ∼0.6 M_⊙. We investigated, for the first time, binaries in the lower part of the MS, with $M_{bin}^1\le 0.5M_{\odot }$, located in the outer regions of NGC 288, one of the five GCs examined by Milone et al. (2020). The outer regions of a GC are particularly interesting because this is where simulations of the dynamical evolution of binaries in multiple-population GCs predict larger differences between the local 1P and 2P binary incidences.

2 Observations and data reduction

We analyzed data collected using the F606W and F814W filters of the Wide Field Channel of the Advanced Camera for Surveys (ACS/WFC), and F110W and F160W filters of the IR channel of the Wide Field Camera 3 (WFC3/IR) on board the Hubble Space Telescope (HST). The ACS/WFC images are part of the program GO 12193 (PI Lee), while the WFC3/IR images are from the program GO 16289 (PI: Milone). Details about these observations are summarized in Table 1.

For data reduction, we performed effective point spread function photometry following the procedure developed by Jay Anderson (Anderson et al. 2008) using the KS2 program. In summary, the process begins by detecting the brightest and most isolated stars in the image. Their fluxes and positions are then determined, after which they are subtracted from the image. Subsequent iterations continue to identify, measure, and subtract progressively fainter stars closer to their neighbors. Three methods are available to achieve optimal photometry and astrometry for different stellar luminosity ranges: method I for relatively bright stars, method II for faint stars, and method III for crowded regions. We adopted method I, which provides flux and position for all stars producing distinct peaks within a 5 × 5 pixel region, using the best Point Spread Function model corresponding to that position. This procedure is repeated for each exposure, and in the end stellar fluxes and positions are properly averaged together. Stellar magnitudes are calibrated to the Vega mag system by computing the aperture correction to the instrumental magnitudes and applying the photometric zero-points available at the space telescope science institute web page (for more details, see Milone et al. 2023). Stellar coordinates are corrected for geometric distortion by using the solutions provided by Anderson et al. (2008) and Anderson (2022). We compared the positions of stars collected at different epochs to derive the proper motions, which we used to identify the probable cluster members and field stars (see, e.g., Anderson & King 2003; Milone et al. 2023, for details). The CMDs m_F814W versus m_F606W − m_F814W and m_F160W versus m_F110W − m_F160W, derived with this procedure, are shown in the left and right panels in Fig. 1, respectively.

To estimate the fraction of 1P and 2P binaries, we compared the observed CMD with simulated CMD constructed using artificial stars (ASs). The ASs test follows the method developed by Anderson et al. (2008). In summary, we generated a sample of ASs equal to ten times the number of observed stars within the instrumental magnitude range of −13.4 to −4.4 in the F814W filter. The F606W, F110W, and F160W magnitudes were calculated using the colors of the fiducial lines of 1P and 2P MS stars. We reduced these ASs using exactly the same procedure as the one adopted for real stars.

3 Binaries along the main sequence

Due to the considerable distance of GCs (e.g., ∼9 kpc for NGC 288, Baumgardt & Vasiliev 2021), the HST is unable to resolve individual components of the binary systems in these stellar clusters. Consequently, binaries appear as single-point sources. The magnitude of an unresolved binary system, m_bin, is determined from the combined fluxes of its components $$m_{\text{bin}}=m_1-2.5\log (1+\frac{F_2}{F_1})$$

with m₁ the magnitude of the primary star (i.e., the most massive one), and F₁ and F₂ the fluxes of the primary and secondary component, respectively.

The procedure used to measure the fraction of binaries with a large mass ratio in NGC 288 is illustrated in the left panel of Fig. 2 and is based on the method by Milone et al. (2012). This method relies on the fact that in the m_F814W versus m_F606W − m_F814W CMD of a GC, binary systems composed of two MS stars are located on the red and bright side of the MS, well separated from the bulk of single MS stars. The exact location of the binaries depends on the mass of the primary star and the mass ratio, q = M₂/M₁, with M₁ and M₂ being the mass of the primary and secondary star, respectively (see Milone et al. 2012; Mohandasan et al. 2024, and references therein).

As illustrated in the left panel of Fig. 2, we identified two regions in the CMD. Region A (gray shaded area in Fig. 2) contains all single MS stars with 19.0 < m_F814W < 23.1 mag, as well as binary systems where the primary star lies within the same luminosity interval. We defined the left boundary (dash-dotted green line) by shifting the MS fiducial line (solid green line) by four times the average color error to the left. The right boundary (dash-dotted red line) corresponds to the equal-mass binaries fiducial (solid red line) shifted by four times the color error to the right. The faint and bright boundaries of the region A (blue lines) correspond to the loci of binary systems with m_F814W = 19.0 and 23.0 mag and q ranging from zero to one. Region B (yellow dashed area in left panel of Fig. 2) is the portion of the region A that includes the binaries with q > 0.5. It is limited on the left side of the CMD by the fiducial of binary systems with q = 0.5 (red dashed line).

The observed stars are shown in the central panel of Fig. 2, where we overlap the lines defined in the left panel onto the m _F814W versus m_F606W − m_F814W CMD. We used yellow starred symbols to mark the probable binaries in the region B of the CMD. The simulated CMD for single MS stars, obtained from ASs, is illustrated in the right panel of Fig. 2.

The fraction of binaries along the MS with q > 0.5 is calculated as $$f_{\mathrm{b}\mathrm{i}\mathrm{n}}^{q>0.5}=\frac{N_{\mathrm{o}\mathrm{b}\mathrm{s}}^B}{N_{\mathrm{o}\mathrm{b}\mathrm{s}}^A}-\frac{N_{\mathrm{s}\mathrm{i}\mathrm{m}}^B}{N_{\mathrm{s}\mathrm{i}\mathrm{m}}^A}$$(1)

where N_obs and N_sim represent the number of observed cluster members and simulated stars, respectively. The resulting binary fraction is $f_{\text{bin}}^{MS}=4.51\pm 0.99\mathrm{\%}$. The uncertainty is derived by means of standard error propagation by associating Poisson errors with the star counts.

Fig. 1 mF814W vs. mF606W − mF814W (left) and mF160W vs. mF110W − mF160W CMD (right) of NGC 288 stars. The error bars to the left of each CMD indicate the typical uncertainties for stars with different magnitudes.

Fig. 2 Procedure for estimating the fraction of binaries among upper MS stars with masses of 0.36−0.75 M⊙. This figure reproduces the mF814W vs. mF606W − mF814W CMD of Fig. 1. The candidate binaries are represented with yellow stars and the remaining stars with small gray circles. The solid green line represents the MS fiducial line, while solid and dashed red lines are the fiducial lines of binaries with mass ratios q = 0.5 and q = 1, respectively. The dashed-dotted green and red lines are the MS fiducial line and the fiducial line of binaries with q = 1 shifted by three times the color error. The blue lines mark the locus of binaries where the mass of the primary component corresponds to 0.36 (upper line) and 0.75 M⊙ (bottom line) and the mass ratio ranges from 0 to 1. In the simulated CMD (right panel), stars located within region B are displayed as gray triangles.

Fig. 3 mF160W vs. mF110W − mF160W CMD (left) and mF606W vs. (mF606W − mF160W) − 0.35(mF110W − mF160W) pseudo-CMD of NGC 288 stars. The probable 1G and 2G M dwarfs, selected from the CMD in the left panel, and the corresponding fiducial lines are colored red and blue, respectively. Dotted lines represent the fiducials of equal-mass binaries (i.e. the MS shifted upward by 0.752 mag). The remaining stars are colored gray, and the fiducial lines of the upper MS are represented by black lines. Typical error bars are shown in red on the left.

Fig. 4 Procedure for selecting and characterizing the binaries used to infer the incidence of 1P and 2P binaries. Left: zoomed-in view of the mF606W vs. (mF606W − mF160W)−0.35(mF110W− mF160W) diagram at the bottom of the MS. Middle: mF160W vs. mF110W−mF160W CMD. Right: verticalized mF160W vs. ΔF110W, F160W diagram (top) and the ΔF110W, F160W cumulative distribution (bottom) for the probable binaries. The probable binaries, which are represented with stars, are selected from the left-panel diagram, where they lie within the yellow-shaded area of the CMD. The solid red and blue lines are the fiducials of 1P and 2P stars, and the dotted lines mark the corresponding loci of equal-mass binaries. See the text for details.

4 Binaries among multiple populations

To estimate the binary fraction among multiple populations, we used the method from Milone et al. (2020), which relies on two photometric diagrams. The first enables the identification of 1P and 2P stars. In the second, the 1P and 2P sequences are nearly superimposed on each other in order to select binaries.

As first photometric diagram, we used the m_F160w versus m_F110w−m_F160W CMD, which has proven effective in distinguishing the multiple populations among M dwarfs (Milone et al. 2014; Dondoglio et al. 2022). This is because the F160W luminosity of M dwarfs is significantly affected by flux absorption due to molecules containing oxygen atoms, whereas the F110W flux is minimally dependent on the stellar oxygen abundance. As a consequence, the 1P stars, which have higher O abundance than 2P stars, exhibit fainter F160W magnitudes and bluer F110W−F160W colors than 2P stars with the same luminosity (e.g., Milone et al. 2012; Dotter et al. 2015).

The presence of multiple populations in NGC 288 is visible in the CMD on the right panel of Fig. 1. Indeed, the color broadening of stars fainter than the MS knee (m_F160W ≳ 20.1 mag) exceeds what can be explained by photometric errors alone, suggesting the presence of stars with different O abundances (see also Dondoglio et al. 2022). We highlighted the probable 1P and 2P stars in red and blue, respectively, in the left panel of Fig. 3, finding that 1P and 2P stars are 59.5 ± 3.8% and 40.5 ± 2.9% of the total M-dwarf population, respectively.

As second photometric diagram, used to select the binaries among multiple populations, we used the m_F606W versus (m_F606W − m_F160W) − 0.35(m_F110W − m_F160W) pseudo-CMD, in which the fiducial lines of 1P and 2P stars overlap (right panel in Fig. 3). In this diagram, we selected a sample of 12 candidate binaries that we used to determine the incidence of 1P and 2P binaries. These stars are located within the dashed yellow area shown in the left panel of Fig. 4 and are represented with starred symbols. The left boundary of this region is defined by the fiducial line of 2P single MS stars, shifted by four times the photometric error to exclude single MS stars with large errors. The right boundary is defined by the fiducial of 2P equal-mass binaries, shifted by four times the color error to account for binaries appearing redder due to photometric errors. We set the upper limit to m_F606W = 22.1 mag because it is not possible to distinguish 1P and 2P stars at brighter magnitudes. The lower limit is set to m_F606W = 23.75 mag because, at fainter magnitudes, it is challenging to disentangle binaries and single stars due to the large observational errors.

In the middle panel of Fig. 4, we display the selected binaries in the m_F160w versus m_F110W−m_F160w CMD, which is used to derive the verticalized diagram m_F160w versus Δ_F110w,F160w. The abscissa of this diagram is calculated in such a way that the two fiducials of equal-mass binaries translate into vertical lines. Specifically, we used the relation $${\mathrm{\Delta }}_{\mathrm{F}110\mathrm{W},\mathrm{F}160\mathrm{W}}=(X-X^{1\mathrm{P}-1\mathrm{P}})/(X^{2\mathrm{P}-2\mathrm{P}}-X^{1\mathrm{P}-1\mathrm{P}})$$(2)

where X is the F110W−F160W color of the selected binary, and X^1P–1P and X^2P–2P are the colors of the fiducial of equal-mass 1P and 2P binaries, respectively, at the same magnitude as the selected binary. The resulting verticalized diagram is shown in the upper right panel of Fig. 4 and is used to derive the cumulative distribution function (CDF) of the selected binaries (lower-right panel of Fig. 4).

To determine the incidence of binaries among 1P and 2P stars, we compared the observations with simulated diagrams that include binaries generated from ASs. The latter exhibit various combinations for the incidence of 1P binaries, 2P binaries, and mixed binary systems, with the fraction of each group of binaries ranging from 0 to 100% in steps of 1%.

We analyzed the simulated diagrams by using the same procedure adopted for the real stars. We derived the CDF of the binaries for each combination of the fraction of the total number of binaries in 1P-1P (f^1P), 2P-2P (f^2P), and 1P-2P (f^mix) mixed binaries and compared it with the corresponding distribution derived for real stars by means of χ². The simulation that provides the best match with the observation corresponds to $f^{1P}={97}_{-3}^{+1}\mathrm{\%},f^{2P}=3\pm 1\mathrm{\%}$, and no mixed binaries $f^{\text{mix}}=0_{-0}^{+1}\mathrm{\%}$.

The uncertainties are derived through a Monte Carlo analysis using 1000 simulated binary samples, each containing the same number of binaries as the observed sample and using the best-fit fractions for 1P, 2P, and mixed binaries. For each simulation, we calculated the values of f^1P, f^2P, and f^mix as we did for the real stars. The upper and lower uncertainties are represented by the 25th and 75th percentiles of the resulting values.

The left and middle panels of Fig. 5 show the m_F606W versus (m_F606W − m_F160W)−0.35(m_F110W − m_F160W) and m_F160W versus m_F110W − m_F160W simulated diagrams, respectively, that provide the best match with the observations. In the right panels we compare the corresponding verticalized m_F160W versus Δ_{F110W, F160W} diagram (top) and the Δ_{F110W, F160W} CDF (bottom) for simulated and observed binaries. The 1P and 2P binary fractions found in our analysis correspond to binary incidences within the 1P(f_{1P, bin}) and 2P(f_{2P, bin}) populations equal to, respectively, f_{1P, bin} = 6.41 ± 1.74% and f_{2P, bin} = 0.29 ± 0.22%.

Fig. 5 Left: simulated diagrams assuming the best-fit values for the fractions of 1P, 2P, and mixed binaries (same layout as Fig. 4). Gray triangles represent artificial binaries within the yellow region, as defined in Section 4. Middle: mF160W vs. mF110W − mF160W CMD, showing the artificial binaries selected in the left-panel diagram. Right: artificial binaries compared with the selected observed binaries, displayed as in Fig. 4.

5 Radial and mass distribution of the binaries in NGC 288

The fact that the fiducial lines of 1P and 2P binaries are nearly coincident in the m_F606W versus (m_F606W − m_F160W) −0.35(m_F110w−m_F160w) plane makes this pseudo-CMD and ideal tool for inferring the overall binary fraction. To do this, we used the method described in Section 2, which is illustrated in Fig. 6, where we compared the observed pseudo-CMD with the simulated pseudo-CMD for single stars. The results are summarized in Table 2 and show that the fraction of binaries with q > 0.5 in the studied field of view is 3.93 ± 0.95%. This value is consistent, within one sigma, with that inferred from the m_F814W versus m_F606W−m_F814W CMD.

Milone et al. (2012) studied the binaries within the half-light radius of NGC 288 and derived a fraction of binaries with q > 0.5 of 5.6 ± 0.4%, which is slightly larger than that observed in our work, at a distance from the center of the cluster of about three half-light radii. The radial behavior of the binary fraction is illustrated in Fig. 7 and is consistent with a flat distribution with some hints of a larger binary fractions in the center.

The m_F606W versus (m_F606W−m_F160w)−0.35(m_F110w− m_F160w) diagram allows us to analyze the binaries with 19.5 < m_F606W < 24 mag, which corresponds to a mass range between 0.39 and 0.75 M_⊙. This mass range is slightly smaller than that studied with the m_F814W versus m_F606W−m_F814W CMD (0.36−0.75 M_⊙). To investigate the relation between the binary fraction and the mass of the primary star, we divided the studied F606W magnitude interval into three bins of 1.5 mag each and calculated the fraction of binaries with q > 0.5 in each bin. As indicated in Table 3, the fractions of binaries in the three luminosity intervals are consistent with each other at one-sigma level. We conclude that there is no evidence of a dependence of the binary fraction on the stellar mass.

Fig. 6 Procedure for estimating the fraction of binaries among lower MS stars with masses of 0.39–0.75 M⊙, following the methodology explained in Section 2. The colors and symbols are consistent with those in Fig. 2.

Fig. 7 Fraction of binaries with q > 0.5 as function of the radial distance from the center of NGC 288 (the radial distance, R, is normalized by the half-light radius, Rhl). Gray circles are measurements from Milone et al. (2012), and black diamond and triangle refer to the binary fraction measured in the magnitude range 19.0 ≤ mF814W ≤ 23 and 19.5 ≤ mF606W ≤ 24, respectively. The dotted blue line marks the core radius (1.35 arcmin). Horizontal and vertical black lines represent the radial coverage of each data point and the uncertainty associated with the measured fraction of binaries, respectively.

6 Simulations

To further investigate the dynamical processes influencing the observed profiles of NGC 288, we conducted a Monte Carlo simulation following the dynamical evolution of a multiple-population GC and its binary population. We would like to emphasize that this simulation is not intended to precisely model NGC 288, but rather to provide insights into the dynamics of multiple-population GCs and their binary stars.

The simulation was performed using the MOCCA code on the Quartz supercomputer at Indiana University. The MOCCA code accounts for the effects of two-body relaxation, binary-single star and single-single star close encounters, stellar evolution, and tidal truncation (see, e.g., Hypki & Giersz 2013; Hypki et al. 2025 for further details about the code). The model begins with 10⁶ stars with stellar masses following a Kroupa (2001) initial mass function between 0.1 M_⊙ and 100 M_⊙ and a 10% binary fraction. The spatial distribution of the 2P and 1P follows King models (King 1966) with central dimensionless potential values, W₀, equal to 7 and 5, respectively. Following the general results of a number of simulations modeling the formation of the 2P (see, e.g., D’Ercole et al. 2008; Bekki 2010; Calura et al. 2019), we start our simulation with the 2P initially centrally concentrated in the inner regions of the 1P system. For the simulation presented here, we adopted an initial ratio of the 2P to the 1P half-mass radius equal to 0.05 (see Vesperini et al. 2021; Livernois et al. 2024 for other simulations starting with this and other values of the 2P concentration; see also Hypki et al. 2022, 2025 for a survey of models exploring various initial conditions). The initial configuration placed 15% of the total initial mass of the cluster in the 2P population (the final fraction of the total mass in 2P stars at 12 Gyr is equal to about 0.51 and similar to the 2P fraction, ∼0.46 ± 0.03, found for NGC288 in Milone et al. 2017b). To align the simulation with the observational constraints, for the results shown in Figures 8 and 9 we restricted the analysis to single stars on the MS with masses between 0.36 and 0.75 M_⊙, and binary stars to MS-MS binaries with mass ratios q > 0.5, where the primary star has a mass between 0.36 and 0.75 M_⊙.

The simulation was evaluated at an evolutionary age of 12 Gyr, consistent with the typical age of old GCs. At 12 Gyr, the fraction of the total mass in 2P stars in our model is about 50%, a value well within the range of those found in Galactic GCs (see, e.g., Milone et al. 2017b).

Each radial profile presented in Figures 8, 9 and 10 was derived from 100 different random realizations of the 2D projection of the 3D data generated by the simulation. For each plot, we display the median profile along with shaded regions representing the 25th and 75th percentiles of the 100 realizations.

One of the results emerging from the observational analysis in the previous sections is the differing degrees of 1P/2P spatial mixing between single stars and binary stars. As shown by Jang et al. (2022) and discussed later in Section 7, the fraction of 2P single stars does not vary with the distance from the cluster’s center, while the fraction of 2P binaries decreases with cluster-centric distance. In the top panel of Fig. 8 we present the radial profiles of the fraction of binaries (solid line) and single stars (dashed line) in the 2P. This figure high-lights a more significant degree of spatial mixing among single stars, reflected in a flatter radial profile, whereas the binary population exhibits a steeper radial gradient starting from smaller radii. This pattern is consistent with binary stars undergoing a slower mixing process compared to single stars. This behavior is consistent with the findings of Hong et al. (2019), who showed that binary stars experience additional dynamical interactions (such as ionization and ejection) that can delay their mixing compared to single stars.

In the lower panel of Fig. 8, we present the radial variation of the ratio between the fraction of binaries and the fraction of single stars in the 2P. A steep radial gradient is observed, further emphasizing that binaries are less mixed than single stars. This result confirms that a GC can exhibit varying degrees of spatial mixing for single and binary stars, consistent with the observational findings for NGC 288. The top panel of Fig. 9 shows the radial profile of binary incidences in each population, while the bottom panel displays the ratio of 1P to 2P binary incidences. In the outer regions, the binary incidence in the 1P population is notably higher than that in the 2P. However, as we move toward the inner regions, the ratio of the incidences tends to 1 (see also the discussion in Milone et al. 2020). The radial variation observed in our simulation aligns with the trend seen in the observational data for NGC 288: Milone et al. (2020) reported similar binary incidences for the 1P and 2P in the cluster’s inner regions, while the results presented in this paper show that the outer regions are characterized by a higher binary incidence in the 1P population compared to the 2P population.

In the top panel of Fig. 9, we also present the global binary incidences for both populations, demonstrating that the 2P global binary incidence is smaller than that of the 1P. This difference arises from the more centrally concentrated spatial distribution and the higher binary ionization rate in the 2P (see Vesperini et al. 2011; Hong et al. 2015, 2016 for the first studies of the dynamics of binaries in multiple-population clusters based on more idealized N-body simulations and Sollima et al. 2022; Hypki et al. 2022, 2025 for more recent and realistic Monte Carlo simulations).

As 1P and 2P single and binary stars mix and interact, binary-binary and single-binary encounters can result in the formation of “mixed binaries”, consisting of one 1P star and one 2P star. The encounter rate, which increases toward the cluster center, plays a key role in the formation of these mixed binaries. Fig. 10 shows that the fraction of the total binary population in mixed binaries in our model is highest near the cluster center and declines steeply in the outer regions (given the small total number of mixed binaries, to better illustrate the radial variation of the mixed binary fraction in this figure we used all the binaries in the model at 12 Gyr and not just those following the observational selection criteria). The radial variation observed in our simulation aligns with the observational trend found in NGC 288 with the data analyzed in this paper and in Milone et al. (2020).

Fig. 8 Top: radial profile of the fraction of binaries (solid line) and single stars (dashed line) belonging to the 2P (the projected distance from the cluster’s center, R, is normalized by the half-light radius, Rhl). Bottom: radial profile of the ratio of the fraction of binaries to the fraction of single stars belonging to the 2P. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection.

Fig. 9 Top: radial profile of the binary incidence in the 1P (red line) and the 2P (blue line). The global binary fractions for each population are indicated by the horizontal lines. Bottom: ratio of the binary incidence in the 1P to that in the 2P. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection.

Fig. 10 Radial profile of the fraction of binaries in mixed binaries. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection. The global fraction of binaries in mixed binaries is indicated by the horizontal dashed line.

7 Summary and discussion

Binary systems provide valuable insight into the formation and dynamical history of multiple stellar populations in GCs. Several formation scenarios suggest that 2P stars form in a high-density subsystem embedded in a less dense 1P cluster; the differences in the 1P and 2P environment set at the time of the cluster formation can affect the survival and evolution of 1P and 2P binary stars. Therefore, the present-day binary fractions of 1P and 2P stars are valuable indicators of the distinct environments in which these stellar populations formed and their subsequent evolution.

Previous photometric studies mainly concentrated on binaries within the innermost regions of GCs, within ∼1 R_hl (Milone et al. 2020), while spectroscopic studies tend to focus on more external regions (D’Orazi et al. 2010; Lucatello et al. 2015; Dalessandro et al. 2018). To fully capture the possible differences in the incidence of 1P and 2P binaries predicted by theoretical studies (see Vesperini et al. 2011; Hong et al. 2015; Sollima et al. 2022; Hypki et al. 2022, 2025), along with the occurrence of mixed binaries (comprising one 1P and one 2P star) formed through interactions with other binaries and single stars, a comprehensive and consistent analysis of the radial variation in 1P and 2P binary frequencies from the cluster center to its outer regions is essential.

The central regions are highly influenced by the cluster’s dynamic evolution; frequent interactions in these dense environments often erase traces of primordial binaries. Conversely, the outer, less dense regions are more effective at preserving this information.

In this study, we examined binaries in the outer region of NGC 288 at ∼5.4 arcmin from the cluster center (∼ 2.5 R_hl) and combined our findings with results obtained by Milone et al. (2020) in the cluster center. To do this, we used HST images collected in the F606W and F814W bands with ACS/WFC, as well as in the F110W and F160W bands with WFC3/IR. To measure the binary fraction among 1P and 2P stars, we adopted the method developed by Milone et al. (2020), thus providing homogeneous investigation of binaries in both the cluster center and in the external field.

First, we disentangled the two stellar populations using the m_F160w versus m_F110w − m_F160w CMD, which effectively disentangles stellar populations with different O abundances below the MS knee. Next, we identified a new photometric diagram where the sequences of the two stellar populations overlap, the m_F606W versus (m_F606W − m_F160W)−0.35(m_F110W − m_F160W) pseudo-CMD, and we selected a sample of 12 binaries. To infer the incidence of 1P and 2P binaries, we compared the CDF of the selected binaries on the verticalized m_F160w versus Δ_F110w,F160W diagram with those derived from a grid of simulated diagrams with different fractions of 1P, 2P, and mixed binaries.

Our analysis reveales a predominance of 1P binaries compared to 2P binaries in the outer regions of NGC 288 where ${97}_{-3}^{+1}\mathrm{\%}$ of the binaries are 1P-1P binaries, 3 ± 1% 2P-2P binaries, and none are mixed binaries. Assuming that 2P stars represent 40.5% of the total stellar population (Fig. 11), we find that the fraction of 1P binaries relative to the total number of 1P stars is approximately 20 times larger than that of 2P stars. Considering a fraction of binaries with q > 0.5 of 3.93 ± 0.95%, we infer that the incidence of binaries is 6.4 ± 1.7% within the 1P and 0.3 ± 0.2% for the 2P.

In the cluster’s center, on the other hand, Milone et al. (2020) measured a similar fraction of 1P and 2P binaries, respectively 46 ± 8% and 40 ± 8% (Fig. 11), and similar binary incidences within each population (with a ratio of the 1P to the 2P binary incidences equal to 1.0 ± 0.3). The radial variation in the 1P and 2P binary incidences is in general agreement with the predictions of models of multiple-population clusters starting with a denser 2P system embedded in the central regions of a less concentrated 1P cluster: the initial differences between the 2P and 1P lead to a more rapid binary disruption for the 2P during the cluster’s dynamical evolution (see Section 6; see also Hong et al. 2015, 2016). Our study also reveals a difference between the degree of spatial mixing of single 1P and 2P stars and that of 1P and 2P binaries; as shown in the upper panel of Fig. 11, the fraction of 2P stars (as measured by Jang et al. 2022) does not vary with the radial distance from the cluster’s center, indicating that the two populations are mixed and lost memory of the initial differences between their spatial distribution. The fraction of the total number of binary stars in 2P binaries, on the other hand, is characterized by a significant radial variation and, as discussed above, decreases in the cluster’s outer regions (bottom panel of Fig. 11). The differences between the degree of mixing of single and binary stars are also in general agreement with the results of the models presented in Section 6, where we have shown that spatial mixing of the two populations occurs on a longer timescale for binaries than for single stars (see Hong et al. 2019 for the first discussion of the various dynamical aspects behind the delayed mixing of 1P and 2P binaries, and Dalessandro et al. 2018 for the possible manifestation of this difference in the kinematics of 1P and 2P stars in NGC 6362). Finally, we point out that our analysis did not find evidence of mixed (1P-2P) binaries in the cluster’s outer regions, while Milone et al. (2020) found a fraction of mixed binaries equal to 14 ± 7% in the central regions. Mixed binaries have been predicted (Hong et al. 2015, 2016) to form in binary-binary and binary-single interactions during which one of the binary components is replaced by a star from a different population. As shown in the simulations presented in Section 6, these binaries are expected to form mainly in the cluster’s inner regions where binary interactions are more frequent and the difference between the inner and outer fraction of mixed binaries is consistent with the predictions of our simulations.

Fig. 11 Top: fraction of 2P stars as a function of the radial distance from the cluster center. Gray circles are measurements from Jang et al. (2022), and black diamond from this research. Horizontal and vertical black lines represent the radial coverage of each data point and the corresponding error bars, respectively. Bottom: fraction of 1P and 2P binaries with respect to the total number of binaries, fbin, as a function of the radial distance from the cluster center. Circle markers are measurements from Milone et al. (2020), and diamond markers from the present work. Red and blue refer to 1P and 2P binary fractions, respectively.

Acknowledgements

This work has received funding from “PRIN 2022 2022MMEB9W – Understanding the formation of globular clusters with their multiple stellar generations” (PI: Anna F. Marino), from INAF Research GTO-Grant Normal RSN2-1.05.12.05.10 – (ref. Anna F. Marino) of the “Bando INAF per il Finanziamento della Ricerca Fondamentale 2022”, and from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 101034319 and from the European Union – NextGenerationEU (beneficiary: T. Ziliotto). E.V. acknowledges support from NSF grant AST-2009193.

References

All Tables

All Figures

	Fig. 1 mF814W vs. mF606W − mF814W (left) and mF160W vs. mF110W − mF160W CMD (right) of NGC 288 stars. The error bars to the left of each CMD indicate the typical uncertainties for stars with different magnitudes.
In the text

Fig. 2 Procedure for estimating the fraction of binaries among upper MS stars with masses of 0.36−0.75 M⊙. This figure reproduces the mF814W vs. mF606W − mF814W CMD of Fig. 1. The candidate binaries are represented with yellow stars and the remaining stars with small gray circles. The solid green line represents the MS fiducial line, while solid and dashed red lines are the fiducial lines of binaries with mass ratios q = 0.5 and q = 1, respectively. The dashed-dotted green and red lines are the MS fiducial line and the fiducial line of binaries with q = 1 shifted by three times the color error. The blue lines mark the locus of binaries where the mass of the primary component corresponds to 0.36 (upper line) and 0.75 M⊙ (bottom line) and the mass ratio ranges from 0 to 1. In the simulated CMD (right panel), stars located within region B are displayed as gray triangles.

In the text

Fig. 3 mF160W vs. mF110W − mF160W CMD (left) and mF606W vs. (mF606W − mF160W) − 0.35(mF110W − mF160W) pseudo-CMD of NGC 288 stars. The probable 1G and 2G M dwarfs, selected from the CMD in the left panel, and the corresponding fiducial lines are colored red and blue, respectively. Dotted lines represent the fiducials of equal-mass binaries (i.e. the MS shifted upward by 0.752 mag). The remaining stars are colored gray, and the fiducial lines of the upper MS are represented by black lines. Typical error bars are shown in red on the left.

In the text

Fig. 4 Procedure for selecting and characterizing the binaries used to infer the incidence of 1P and 2P binaries. Left: zoomed-in view of the mF606W vs. (mF606W − mF160W)−0.35(mF110W− mF160W) diagram at the bottom of the MS. Middle: mF160W vs. mF110W−mF160W CMD. Right: verticalized mF160W vs. ΔF110W, F160W diagram (top) and the ΔF110W, F160W cumulative distribution (bottom) for the probable binaries. The probable binaries, which are represented with stars, are selected from the left-panel diagram, where they lie within the yellow-shaded area of the CMD. The solid red and blue lines are the fiducials of 1P and 2P stars, and the dotted lines mark the corresponding loci of equal-mass binaries. See the text for details.

In the text

Fig. 5 Left: simulated diagrams assuming the best-fit values for the fractions of 1P, 2P, and mixed binaries (same layout as Fig. 4). Gray triangles represent artificial binaries within the yellow region, as defined in Section 4. Middle: mF160W vs. mF110W − mF160W CMD, showing the artificial binaries selected in the left-panel diagram. Right: artificial binaries compared with the selected observed binaries, displayed as in Fig. 4.

In the text

	Fig. 6 Procedure for estimating the fraction of binaries among lower MS stars with masses of 0.39–0.75 M⊙, following the methodology explained in Section 2. The colors and symbols are consistent with those in Fig. 2.
In the text

Fig. 7 Fraction of binaries with q > 0.5 as function of the radial distance from the center of NGC 288 (the radial distance, R, is normalized by the half-light radius, Rhl). Gray circles are measurements from Milone et al. (2012), and black diamond and triangle refer to the binary fraction measured in the magnitude range 19.0 ≤ mF814W ≤ 23 and 19.5 ≤ mF606W ≤ 24, respectively. The dotted blue line marks the core radius (1.35 arcmin). Horizontal and vertical black lines represent the radial coverage of each data point and the uncertainty associated with the measured fraction of binaries, respectively.

In the text

Fig. 8 Top: radial profile of the fraction of binaries (solid line) and single stars (dashed line) belonging to the 2P (the projected distance from the cluster’s center, R, is normalized by the half-light radius, Rhl). Bottom: radial profile of the ratio of the fraction of binaries to the fraction of single stars belonging to the 2P. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection.

In the text

	Fig. 9 Top: radial profile of the binary incidence in the 1P (red line) and the 2P (blue line). The global binary fractions for each population are indicated by the horizontal lines. Bottom: ratio of the binary incidence in the 1P to that in the 2P. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection.
In the text

	Fig. 10 Radial profile of the fraction of binaries in mixed binaries. The shaded regions represent the 25th and 75th percentiles of the 100 realizations of the 2D random projection. The global fraction of binaries in mixed binaries is indicated by the horizontal dashed line.
In the text

Fig. 11 Top: fraction of 2P stars as a function of the radial distance from the cluster center. Gray circles are measurements from Jang et al. (2022), and black diamond from this research. Horizontal and vertical black lines represent the radial coverage of each data point and the corresponding error bars, respectively. Bottom: fraction of 1P and 2P binaries with respect to the total number of binaries, fbin, as a function of the radial distance from the cluster center. Circle markers are measurements from Milone et al. (2020), and diamond markers from the present work. Red and blue refer to 1P and 2P binary fractions, respectively.

In the text