© The Authors 2026

1 Introduction

Open clusters (OCs) serve as fundamental astrophysical laboratories for testing theories of stellar evolution and stellar dynamics. As OCs age, they undergo continuous dynamical evolution driven by two-body relaxation, leading to energy equipartition and mass segregation (Spitzer 1987). In dynamically mature clusters, massive stellar components are expected to efficiently sink toward the central core. Consequently, the presence of massive stars in the extended cluster halo presents a compelling paradox, hinting at violent dynamical mechanisms capable of overcoming the deep central potential well.

Blue straggler stars (BSSs), first identified by Sandage (1953), are enigmatic stellar products that appear younger, hotter, and more massive than the main-sequence turnoff (MSTO) of their host cluster. Since standard single-star evolution cannot explain their necessary mass gain, two primary formation channels have been widely accepted (Rain et al. 2024): (i) stable or unstable mass transfer (MT) or coalescence within primordial binary systems (McCrea 1964), and (ii) direct stellar collisions mediated by few-body dynamical interactions (Hills & Day 1976). Recent high-precision photometric studies have further revealed that BSSs often form distinct sequences in the color-magnitude diagram (CMD) – a “blue” sequence primarily associated with collisional products, and a “red” sequence predominantly formed via MT (Milone et al. 2012). While collisions are prevalent in the high-density environments of globular clusters, the relatively sparse environments of OCs have traditionally been viewed as laboratories dominated by the MT channel (Leigh et al. 2011). Recent work by Ferraro et al. (2026) further emphasizes that the efficiency of binary-related formation channels is strongly modulated by the host cluster’s environmental conditions, such as density and age, providing a unified framework for understanding the diverse BSS populations observed across different stellar systems.

NGC 188 (C 0039+850 in the IAU nomenclature), with an age of approximately 7 Gyr (Cantat-Gaudin et al. 2020; Vaidya et al. 2020; Penev & Schussler 2022; Alvarez-Baena et al. 2024; Sun et al. 2025; Yakut et al. 2025; Narayan et al. 2026), is one of the oldest and most well-studied OCs in the Milky Way, serving as a fundamental benchmark for BSS evolution. Extensive radial velocity surveys by the WIYN Open Cluster Study (WOCS) have revealed that a vast majority of its central BSS population resides in long-period binaries, pointing toward a formation history heavily dominated by MT (Mathieu & Geller 2009; Geller & Mathieu 2011; Sun et al. 2022). Consequently, NGC 188 has often been cited as a “quiescent” cluster where violent dynamical encounters play a negligible role.

However, modern N-body simulations suggest a bifurcated reality: even in relatively low-density clusters, the core can act as a persistent “dynamical furnace” (Hurley et al. 2005). Close encounters – such as binary-single or binary-binary scattering – can trigger BSS formation via collisions or interaction-enhanced MT. A critical yet historically elusive signature of such encounters is the kinetic recoil. The conservation of momentum during a strong gravitational interaction dictates that the reaction products will receive a significant “kick” velocity (Perets 2009; Purohit et al. 2024). If this kick exceeds the local velocity dispersion but remains below the cluster escape velocity, the BSS will be ejected into the cluster halo on a highly eccentric orbit.

Despite these theoretical predictions, observational evidence of such “escaped” BSSs in OCs has been exceedingly scarce. Most past surveys have been spatially biased toward the high-density cluster core, leaving the expansive halo largely unexplored (Geller et al. 2008; Cantat-Gaudin et al. 2018). The advent of Gaia Data Release 3 (DR3) (Gaia Collaboration 2023) now provides the unprecedented astrometric precision required to conduct a “kinematic forensic” analysis of cluster outskirts, searching for high-velocity members that retain the dynamical memory of their violent birth.

In this work, we present a comprehensive census of the stellar population in the old OC NGC 188, utilizing the high-precision astrometry provided by Gaia DR3. Our primary objective is to construct a spatially extended membership catalog that reliably traces the cluster structure from the core out to the distant halo. By characterizing the blue straggler population within this extended framework, we aim to constrain their formation mechanisms and specifically investigate the kinematic signatures of dynamical interactions. We focus on identifying high-velocity candidates in the cluster periphery to test whether dynamical ejection processes remain active in evolved cluster environments.

The paper is organized as follows. Section 2 details the data retrieval and the probabilistic membership determination using the UPMASK algorithm. In Section 3, we validate our catalog and assess its completeness relative to recent literature. Section 4 describes the differential reddening correction applied to refine the CMD. In Section 5, we derive the structural parameters and analyze the global dynamical state of the cluster, including mass segregation. The identification and photometric characterization of the blue straggler population are presented in Section 6. Section 7 focuses on the detailed kinematic analysis, including the alignment test and orbit reconstruction for high-velocity candidates. Finally, we summarize our conclusions in Section 8.

2 Data and membership selection

2.1 Gaia DR3 data retrieval

We initiated our study by retrieving all sources within a 1.5° radius centered on the coordinates of NGC 188 (α = 11.817°, δ = +85.255°) (reported in Nizovkina et al. 2025) from the Gaia DR3 (Gaia Collaboration 2023). To ensure the reliability of the kinematic analysis, we applied a series of stringent quality filters: (i) a parallax signal-to-noise ratio ϖ/σ_ϖ ≥ 3 to eliminate poorly constrained distances; (ii) proper motion uncertainties proper motion uncertainties μ_α* < 1 and μ_δ < 1 mas yr⁻¹ mas yr⁻¹ to maintain high-precision velocity vectors; and (iii) a physical parallax constraint ϖ > 0. After these quality cuts, a clean sample was prepared for membership assignment and yielded 9542 sources.

2.2 The modified UPMASK algorithm

To identify cluster members without relying on theoretical isochrones (thereby avoiding circular reasoning for BSS identification), we employed a modified version of the Unsupervised Photometric Membership Assignment in Stellar Clusters (UPMASK) algorithm (Krone-Martins & Moitinho 2014). This unsupervised, data-driven approach determines membership probabilities based on the clustering of stars in the 3D astrometric subspace (${\mu }_{\alpha }^{\ast }$, μ_δ, ϖ).

The core of our implementation consists of an iterative Monte Carlo process. In each of the N_iter = 100 iterations, we perturbed the astrometric parameters of each star within their observational error bars: $$X_{i,\text{perturbed}}=X_{i,\text{obs}}+\mathcal{N}(0,{\sigma }_{X_i}),$$(1)

where X ∈ {${\mu }_{\alpha }^{\ast }$, μ_{δ, ϖ}} represents the astrometric observables, and 𝒩(0, σ_Xi) corresponds to the normal distribution centered on 0 with a dispersion equal to the measurement uncertainty σ_Xi. To ensure that all dimensions contribute equally to the distance-based partitioning, the perturbed observables (ϖ, ${\mu }_{\alpha }^{\ast }$, μ_δ) were first standardized using a StandardScaler, which rescales the data to have zero mean and unit variance. These standardized data were then subjected to K-means clustering, where we partitioned the sample into k = 15 distinct groups. By selecting a relatively high number of groups (k = 15), we over-partition the feature space; this facilitates the isolation of the physical stellar cluster (NGC 188) as a single, high-density entity, preventing it from being numerically merged with the diffuse and more populous background of field stars in the velocity–parallax space. The initial priors for the cluster center in proper-motion space were set to (${\mu }_{\alpha }^{\ast }$, μ_δ) = (−2.3, −0.9) mas yr⁻¹.

These values were adopted based on the established mean kinematics of NGC 188 reported in recent Gaia-based catalogs (e.g., Cantat-Gaudin et al. 2020). Utilizing these informed priors ensures that the membership algorithm (UPMASK) efficiently isolates the cluster population from the field stars by centering the search on the known velocity of the system. For each iteration, the cluster whose centroid was closest to the prior cluster mean (${\overline{\mu }}_{\alpha }^{\ast }\approx -2.3,{\overline{\mu }}_{\delta }\approx -0.9$ mas yr⁻¹) was identified as the “member group.”

A star’s membership probability, P_mem, is defined as the frequency with which it is assigned to the member group over all iterations: $$P_i=\frac{1}{k}\sum _{j=1}^k{I}_{i,j},$$(2)

where I_i,j is an indicator function that takes the value 1 if the star i is assigned to a high-density group (representing the physical cluster) in the j-th iteration of the K-means clustering, and 0 otherwise. Here, k represents the total number of iterations performed on the perturbed datasets. This frequentist approach allows us to transform the discrete clustering outputs into a continuous membership probability, P_i. Our analysis identified 1652 stars with P_mem > 0.5. From this sample, a highly purified subset of 1214 “high-confidence” members (P_mem > 0.9) was used to define the cluster’s fiducial sequences in the CMD. This robust membership list provides the foundation for identifying BSSs that have been dynamically scattered into the cluster halo.

3 Comparison with previous catalogs

To validate the reliability of our membership selection, we cross-matched our results with the recently published OC catalog by Hunt & Reffert (2024) (hereafter Hunt24). Figure 1 illustrates the comparative analysis between our sample (N = 1652) and the Hunt24 catalog (N = 1171).

The vector voint viagram (VPD) demonstrates that both samples are centered on the same kinematic mean, yet our UPMASK-derived sample captures a broader range of the cluster’s velocity dispersion in the outer regions. The most striking difference is observed in the spatial distribution (second panel of Figure 1); whereas Hunt24 is largely restricted to the central 0.25°, our pipeline successfully identifies members out to the 0.5° tidal radius limit. The best fitting Padova isochrone (dotted black line) is given in right panel of Figure 1.

The parallax distributions show excellent agreement between the two catalogs. Taking advantage of our high-purity member sample, we derive a precise macroscopic distance to NGC 188. Based on the N members identified in this work, we obtain a mean parallax of ϖ = 0.5462 ± 0.0035 (stat.) ± 0.015 (sys.) mas, where the systematic uncertainty accounts for the Gaia astrometric zero-point calibration floor for localized fields (Lindegren et al. 2021). Simple inversion of this high-precision parallax yields a cluster distance of d = 1830.9 ± 11.6 (stat.) ± 50.3 (sys.) pc. Photometrically, the additional members identified in this work (blue points in Figure 1) are indistinguishable from the core members in the CMD, following the 7 Gyr isochrone with minimal scatter. This high-purity, high-completeness sample, alongside the newly derived cluster parameters, is critical for the subsequent kinematic forensic analysis of the BSS population.

Fig. 1 Four-panel comparison between our membership results (blue, N = 1652) and the Hunt & Reffert (2024) catalog (orange, N = 1171). Left: vector Point Diagram in proper motion space. Second: spatial distribution in RA and Dec. Third: histogram of Gaia DR3 parallaxes. Right: G vs. (GBP − GRP) CMD. Our method identifies a significant population of halo members that were previously overlooked.

Fig. 2 Differential reddening correction diagnostics. Left: raw CMD with the fiducial MS ridge line (red). Middle: spatial map of the color residuals δE(BP − RP), indicating a fairly uniform extinction field with subtle gradients. Right: final corrected CMD showing sharpened stellar sequences.

4 Differential reddening correction

To maximize the photometric resolution of the CMD, we applied a differential reddening correction based on the multi-step interpolation technique of Milone et al. (2012). Given the 7 Gyr age of NGC 188, a sharp definition of the main sequence (MS) and Blue Straggler region is paramount for identifying low-mass members and subtle kinematic outliers.

We defined a reference sample of “pure” MS stars by selecting high-probability members (P_mem > 0.8) within the range 15.0 < G < 19.0. To prevent unresolved binaries from skewing the fiducial line, we utilized a robust binning strategy, retaining only stars within 2.5σ of the median color in each 0.4 mag interval. The fiducial ridge line was then modeled using a univariate spline (s=0.005).

The local extinction residual for each star, δE(BP − RP), was estimated by adopting the median color residual of its k = 35 nearest neighbors. As illustrated in Figure 2, the spatial distribution of extinction is relatively flat, with minor fluctuations across the cluster halo. We adopted a Gaia DR3 extinction coefficient of A_G/E(BP − RP) = 2.05 to correct both magnitudes and colors.

The success of the procedure is evidenced by the narrowing of the MS sequence: the standard deviation of the color distribution improved from σ_raw = 0.0429 mag to σ_corr = 0.0424 mag. This refinement, though small, ensures that the BSS candidates analyzed in Section 6 are not photometric artifacts resulting from patchy foreground dust.

5 Cluster structure and dynamical state

5.1 Surface density profile and structural parameters

To characterize the spatial structure of NGC 188, we derived the stellar surface density profile using the cleaned membership catalog. As shown in Figure 3 (left), we fit the observed radial distribution with a semiempirical King model (King 1962): $$\rho (r)={\rho }_0{[\frac{1}{\sqrt{1+{\left(r/R_c\right)}^2}}-\frac{1}{\sqrt{1+{\left(R_t/R_c\right)}^2}}]}^2+{\rho }_{\mathrm{b}\mathrm{g}}.$$(3)

In this expression, ρ(r) represents the projected stellar surface density at a radial distance, r, from the cluster center. The parameters ρ₀ and ρ_bg denote the central surface density and the constant background density, respectively. R_c signifies the core radius, while R_t represents the tidal radius, the distance at which the cluster’s density is truncated by the Galactic tidal field.

Our best-fit parameters yield a core radius of R_c = 2.00 ± 0.15 pc and a tidal radius of R_t = 51.2 ± 4.3 pc. The concentration parameter c = log₁₀(R_t/R_c) ≈ 1.41 suggests that NGC 188 is a moderately concentrated system. The smooth fit to the surface density out to the halo regions (r > 15 pc) provides a robust spatial baseline for identifying ejected blue straggler candidates.

Fig. 3 Structural and dynamical properties of NGC 188. Left: surface density profile derived from the GAIA DR3 membership catalog. The solid line represents the best-fit King model, yielding a core radius, Rc = 2.00 ± 0.15 pc, and a tidal radius, Rt = 51.2 ± 4.3 pc. The background density level, ρbg, is subtracted to highlight the cluster’s extent into the halo. Right: mass segregation ratio, ΛMSR, as a function of the number of the most massive stars, NMST. The significant deviation from unity (ΛMSR > 1.2 for the top 50 massive stars) confirms a strong mass segregation effect, indicating that the cluster is dynamically well relaxed because of its 7 Gyr age.

5.2 Structural parameters and potential scale

To establish a self-consistent dynamical model for NGC 188, we reconciled the observed stellar mass with the cluster’s global dynamical properties. We estimated the mass of each member star using a Monte Carlo approach following the methodology of Almeida et al. (2023) and Chi et al. (2026). Based on the fundamental parameters of NGC 188, we generated synthetic stellar populations. For each observed star, we identified its closest counterpart in a synthetic cluster within multiband magnitude space. A star was classified as a binary if its best-matching synthetic counterpart was flagged as a binary, and was assigned the total system mass. This procedure was repeated over 10 000 synthetic realizations, with the final mass (and companion mass, if applicable) taken as the ensemble mean and the uncertainty as the standard deviation. The raw photometric mass derived from our Gaia DR3 membership catalog, M_obs ≈ 1425 M_⊙, represents a lower limit to the true population due to the intrinsic completeness limits of the survey (G ~ 20.5) and the presence of unresolved binaries. Through our integrated mass analysis based on method of Chi & Wang (2025), we derive a total cluster mass for NGC 188 of M_tot ≈ 2180 M_⊙. Following the dynamical mass estimates of Section 5.4, we adopted a total cluster mass of M_tot = 2180 M_⊙ for our gravitational potential, implying a global correction factor of f_corr ≈ 1.53. This factor accounts for both the missing low-mass M-dwarfs and the hidden mass in binary systems, which are particularly prevalent in this cluster (Geller & Mathieu 2012).

By analyzing the probability-weighted cumulative mass profile M(< R), we derived a projected half-mass radius of r_h = 5.59 pc (see Figure 4). This value is slightly more extended than previous estimates restricted to the core (Bonatto et al. 2005; Wang et al. 2015; Gao & Fang 2022), reflecting the high efficiency of our UPMASK-based approach in recovering halo members. For the subsequent orbital back-tracing, we adopted a Plummer potential: $$\mathrm{\Phi }(r)=-\frac{G{M}_{\mathrm{c}\mathrm{l}}}{\sqrt{r^2+a^2}},$$(4)

where the scale parameter a is set by the projected half-mass radius (a = r_h = 5.59 pc), ensuring the mass distribution in our model matches the observed stellar density profile. The structural consistency between our derived r_h and historical benchmarks validates the integrity of our expanded membership catalog and provides a robust dynamical foundation for interpreting the runaway blue straggler population.

5.3 Mass segregation and dynamical maturity

Given the ancient age of NGC 188 (~7 Gyr), significant dynamical relaxation is expected. To quantify the degree of mass segregation, we employed the minimum spanning tree (MST) method (Allison et al. 2009). An MST is defined as the unique set of edges that connects a given set of points (stars) such that the total edge length is minimized without forming any closed loops. Unlike traditional radial density profiles, the MST method is advantageous as it does not require the prior determination of the cluster center and is robust against outliers in the stellar distribution.

The degree of mass segregation is quantified by the mass segregation ratio, Λ_MSR, which is defined as $${\mathrm{\Lambda }}_{MSR}(N)=\frac{⟨{\ell }_{random}⟩}{{\ell }_{massive}}\pm \frac{{\sigma }_{random}}{{\ell }_{massive}},$$(5)

where ℓ_massive is the total edge length of the MST connecting the N most massive stars in the sample. The parameter ⟨ℓ_random⟩ is the average edge length of the MSTs constructed from 500 sets of N randomly selected stars from the entire cluster population, and σ_random is the associated standard deviation.

The physical interpretation of Λ_MSR is straightforward: if Λ_MSR ≈ 1, the spatial distribution of massive stars is indistinguishable from that of the average population (no mass segregation). If Λ_MSR > 1, the massive stars are significantly more concentrated than a random sample of stars, indicating the presence of mass segregation. As presented in Figure 3 (right panel), the Λ_MSR was calculated for the N_MST most massive stars (primarily BSSs and red giants) relative to a random reference population, providing a clear signature of the dynamical state of NGC 188.

We find a significant mass segregation effect, with Λ_MSR ≈ 1.3–1.5 for the top 50 most massive members. This elevation above unity (Λ_MSR > 1) indicates that the most massive stars are more centrally concentrated than the average cluster members. This result confirms that NGC 188 is dynamically “relaxed” or “mature.” The concentration of these massive BSSs in the core increases the probability of multi-body interactions (e.g., 2 + 2 or 1 + 2 scattering), providing a consistent dynamical environment for the ejection of high-velocity BSSs into the cluster halo, as discussed in Section 7.

Fig. 4 Probability-weighted cumulative mass profile for NGC 188. The solid dark curve represents the cumulative stellar mass M(< R) derived from our corrected catalog (Mtot = 2180 M⊙). The intersection of the dashed red lines identifies the projected half-mass radius rh = 5.59 pc. The shaded region denotes the volume containing 50% of the cluster’s total dynamical mass, providing the structural basis for the adopted Plummer potential.

5.4 Mass distribution and relaxation state

Based on the estimated total cluster mass of M_tot ≈ 2180 M_⊙ (Section 5.2) for NGC 188, and taking a characteristic half-mass radius of r_h = 5.59 pc, the calculated half-mass relaxation time is t_rh ≈ 125.4 Myr. Given the cluster’s canonical age of 7 Gyr (logt = 9.85), NGC 188 has survived for approximately 50.6 relaxation times (T_age/t_rh ≈ 50.6).

This ratio identifies NGC 188 as an extremely dynamically mature system, significantly exceeding the threshold typically required for complete mass segregation (T_age/t_rh > 10). In such a state, the cluster’s massive members, particularly binary systems and BSSs, have efficiently subsided into the high-density core. As shown in the CMD (Figure 5), the BSS population is clearly identified above the MSTO, consistent with a history of mass transfer or stellar mergers.

Fig. 5 Color–magnitude diagram of NGC 188 color-coded by estimated stellar mass. The red asterisks mark the identified runaway BSS candidates (see Section 7).

6 Characterization of the blue straggler population

The identification and classification of BSSs in NGC 188 were performed using a rigorous selection procedure based on both theoretical models and empirical cluster sequences. To define the BSS search region, we first fit the cluster’s clean MS and turnoff (TO) with a theoretical PARSEC isochrone (Bressan et al. 2012; Chi et al. 2025) assuming an age of 7.0 Gyr, solar metallicity ([Fe/H] = 0.0), and a distance modulus adjusted to our derived distance.

The BSS selection box was constructed in the de-reddened G₀ versus (G_BP − G_RP)₀ CMD using the following criteria:

1. Membership constraint: Only stars with a Gaia-based membership probability P > 0.9 were considered to minimize field star contamination;

2. Luminosity and color boundaries: Candidates were required to be brighter than the MSTO (G₀ < 15.7) and bluer than the intrinsic color of the TO point ((G_BP − G_RP)₀ ≲ 0.78 mag);

3. Empirical reference sequences: To prevent contamination from the photometric binary sequence and the subgiant branch, we defined a red boundary for the BSS box by shifting the 7.0 Gyr isochrone by −0.05 mag in color. The blue boundary was defined by the zero age main sequence (ZAMS), ensuring that candidates remain within physically plausible limits for hydrogen-burning stars.

Following these criteria, we identified 45 BSS candidates (Figure 6a). This sample includes “classical” BSSs located well above the TO. The resulting sample is consistent with previous catalogs (e.g., Geller et al. 2008; Gosnell et al. 2015) but benefits from the improved astrometric precision and completeness of our Gaia-based membership analysis. Our identified sample of 45 BSS candidates is larger than the 22 objects reported for NGC 188 by Rain et al. (2021). This discrepancy primarily arises from the differences in membership selection and CMD boundary definitions. While Rain et al. (2021) utilized a more conservative, automated selection algorithm based on Gaia DR2 to ensure a high-purity sample across a large number of clusters, our study benefits from the improved astrometric precision of Gaia DR3 and a more inclusive selection box tailored specifically to the multiwavelength properties of NGC 188.

Fig. 6 Comprehensive analysis of the BSS population. Panel a: CMD showing the BSS selection box relative to the Padova isochrone and the empirical M67 ZAMS. Panel b: histogram of color offsets from the ZAMS, defining the blue and red sequences. Panel c: multiparameter binary identification matrix. Circle sizes are proportional to square of RUWE, and colors represent the image morphology (ipd_frac_multi_peak).

6.1 Sequence bifurcation

To distinguish between different formation histories, we compared the BSS distribution against an empirical ZAMS derived from M67, appropriately shifted for the distance and reddening of NGC 188. We define the color residual Δ(G_BP − G_RP)₀ relative to this ZAMS. The population exhibits a bifurcated structure. using a threshold of Δ = 0.08 mag, corresponding to the local minimum in the color residual distribution (Figure 6b), we classify the sample into a Blue Sequence (likely collision products or rejuvenated single stars) and a Red Sequence (likely mass-transfer binaries).

6.2 Binary content and Gaia diagnostics

We utilized Gaia DR3 astrometric and morphological flags to probe the binary nature of the BSS sample. Specifically, we employed the Renormalized Unit Weight Error (RUWE) as a proxy for astrometric binaries and the ipd_frac_multi_peak parameter to identify unresolved companions with significant flux contributions.

Figure 6c displays the results of this multiparameter diagnostic. While a subset of BSSs shows high RUWE values (>1.4) indicative of significant orbital motion, a substantial fraction remains astrometrically “quiet.” This observation supports previous findings that many BSS companions in NGC 188 are low-mass white dwarfs, which exert minimal astrometric pull on the primary BSS (Geller & Mathieu 2011).

6.3 Origin of the bifurcated sequences: Mass transfer versus collisions

The blue straggler population in NGC 188 exhibits a clear preference for the red sequence, which accounts for 88.9% (40 out of 45) of the total identified members. This overwhelming dominance provides a key diagnostic of the cluster’s internal environment.

Formation channels: in an ancient OC such as NGC 188 (~7 Gyr), the stellar density in the core is significantly lower than that of massive globular clusters. Consequently, the rate of direct stellar collisions – the primary architect of the blue sequence – is suppressed. The high fraction of red sequence BSSs suggests that stable mass transfer in primordial binary systems is the primary formation channel. This is consistent with the long-term secular evolution of NGC 188, where binary systems have had gigayears to initiate Roche-lobe overflow as the primary star evolves off the MS.

Spatial distribution and statistical limitations: We investigated whether these two formation channels lead to distinct spatial footprints. Collision-induced BSSs (blue sequence) are theoretically expected to be more massive and thus more centrally concentrated due to dynamical friction.

While the BSSs are predominantly found in the inner regions, our quantitative analysis using a two-sample Kolmogorov–Smirnov (K–S) test yielded a statistic of D = 0.250 with a p value of 0.902. This result indicates that we cannot statistically reject the null hypothesis that both sequences share the same radial distribution. However, we note that this lack of statistical significance is primarily driven by the extremely small sample size of the blue sequence (N = 5). As the cluster continues to evaporate and undergo mass segregation, the dynamical distinction between these populations may become more pronounced, but currently, the mass transfer channel remains the ubiquitous source of BSSs across the entire cluster extent.

7 Kinematic forensics and ejection scenarios

To differentiate between in situ halo BSSs and those dynamically ejected from the core, we performed a rigorous kinematic analysis. We established a kinematic reference frame based on the median proper motion of the cluster’s red giant branch (RGB) stars. The use of RGB stars as the reference population is a standard and robust technique for several key reasons:

1. Astrometric precision: As one of the brightest stellar populations in the cluster, RGB stars have high signal-to-noise measurements in Gaia, leading to the smallest uncertainties in their proper motions (${\mu }_{\alpha }^{\ast }$, μ_δ);

2. High purity sample: The RGB forms a distinct and narrow sequence in the CMD, allowing for an unambiguous selection of member stars with minimal contamination from field objects;

3. Dynamical stability: Being a massive and evolved population, RGB stars are expected to be dynamically relaxed and serve as excellent tracers of the cluster’s gravitational potential well. They constitute a kinematically “cold” population with a low velocity dispersion, providing a stable benchmark against which the motions of other, potentially more dynamic, populations can be measured.

The median motion of this clean RGB sample yields the cluster’s bulk velocity (${\mu }_{\alpha }^{\ast }$, μ_δ = −2.31, −0.92 mas yr⁻¹), a reference frame that is robust against outliers and ideal for identifying stars with peculiar kinematics, such as ejected BSSs. Having characterized the overall BSS population, we now focus on a subset of kinematic outliers located in the cluster halo. Among the 45 candidates, two specific objects (hereafter BSS-E1 and BSS-E2) stand out due to their anomalous velocities and trajectories, suggesting a nonstandard dynamical history. While the full 3D orbital trajectory cannot be uniquely determined without precise radial velocities (RVs), our 2D kinematic analysis on the tangent plane remains robust and conclusive for two fundamental reasons:

1. Lower limit on escape energy: The relative transverse velocities (v_rel) of our BSS candidates, derived strictly from Gaia proper motions, are 5.9 and 8.9 km s⁻¹. In the potential well of NGC 188, the local escape velocity is estimated to be ~2–3 km s⁻¹. Since the total 3D velocity is defined as $v_{\text{tot}}=\sqrt{v_{\text{rel}}^2+v_R^2}$, the inclusion of any nonzero radial component v_R would only strictly increase the total velocity, further cementing their status as unbound, escaping objects;

2. Improbability of chance alignment: The alignment score (cos θ > 0.975, see Section 7.1) provides a powerful geometric constraint. Given the cluster’s high declination (δ ≈ 85°), we have rigorously applied spherical offsets (Δα cos δ) to account for projection effects. The probability that a random field interloper, or even an unperturbed halo member, would possess a velocity vector that precisely back-traces to the sub-parsec cluster core over a ~20 arcmin baseline is vanishingly small.

Therefore, the combination of “supersonic” transverse motion and extreme radial pointing serves as a definitive “kinematic fingerprint” of a recent core-driven ejection event, independent of the unknown radial component.

Fig. 7 Kinematic map of the BSS population in NGC 188. Arrows represent relative proper motion vectors scaled to 8 times their annual displacement. The crimson stars denote the two identified ejected candidates, characterized by high velocities and high alignment scores (cos θ). The core and halo limits are indicated by dashed and dotted circles, respectively. Black crosses denote excluded interlopers.

7.1 Alignment analysis

The key metric used to identify ejection events is the alignment score, cos θ. We define θ as the angle between the star’s position vector (r), measured from the cluster center, and its relative proper motion vector (Δμ = μ_star − μ_cluster). Mathematically, this is expressed as $$\cos \theta =\frac{\mathit{r}\cdot \mathrm{\Delta }\mathit{\mu }}{|\mathit{r}||\mathrm{\Delta }\mathit{\mu }|},$$(6)

where μ_cluster is the median motion of the cluster derived from the RGB members. A value of cos θ ≈ 1 indicates that the star is moving almost perfectly radially outward from the core.

The relative tangential velocities, v_rel, were derived from Gaia DR3 proper motions and the cluster distance. We first calculated the relative proper motion magnitude $\mathrm{\Delta }\mu =\sqrt{{({\mu }_{\alpha }^{\ast }-{\overline{\mu }}_{\alpha }^{\ast })}^2+{({\mu }_{\delta }-{\overline{\mu }}_{\delta })}^2}$. This was then converted into a physical velocity (in kilometers per second) using the relation $$v_{\text{rel\hspace{0.17em}}}=4.7405\times d\times \mathrm{\Delta }\mu ,$$(7)

where d is the distance to NGC 188, and 4.7405 is the unit conversion factor.

In Figure 7, we applied a more stringent kinematic pruning to this high-confidence sample. The “interlopers” (marked as stars with v_rel > 10 km s⁻¹) are objects that passed the initial astrometric cut but are dynamically inconsistent with being bound members of NGC 188. Given the cluster’s low internal velocity dispersion (σ_v ≈ 0.5 km s⁻¹), a relative velocity of 10 km s⁻¹ represents a 20σ deviation. These stars are likely field stars whose proper motions and parallaxes coincidentally overlap with the cluster’s mean values within the Gaia uncertainties.

Specific cases mentioned:

• The star at (−25, +5): This object exhibits a large velocity vector but an alignment score of cos θ ≈ 0.42. Its presence in the halo, combined with its nonradial “random” trajectory, identifies it as a kinematic interloper (field contaminant) rather than a candidate for dynamical ejection.

• High-velocity stars between the dashed circles: While several BSS candidates in the halo region (R_h < r < 2R_h) show relatively large velocities, they were not classified as “Ejected” because they lack the requisite radial alignment (cos θ < 0.9). In a cluster environment, high-velocity stars that do not point away from the core are more likely to be interlopers or stars in the process of being tidally stripped by the Galactic potential, rather than products of core dynamical encounters.

By explicitly plotting these interlopers, we demonstrate the necessity of using both velocity magnitude and the alignment score (cos θ) to isolate true dynamical ejection events from the background noise. As shown in Figure 7 and Table 1, while most halo BSSs (gray points) exhibit random orbital phases, two stars (red stars) display extreme radial trajectories. By “extreme”, we refer to a kinematic state that is statistically and dynamically inconsistent with the relaxed cluster population:

• Directional extremity: Their alignment scores (cos θ ≥ 0.98) imply a nearly “pure” radial motion, which is highly improbable for stars in standard stochastic orbits but characteristic of a “slingshot” ejection from the core;

• Velocity extremity: Their velocities (5.9–8.9 km s⁻¹) are significantly higher than the typical 1D velocity dispersion of NGC 188 (σ_v ≈ 0.5 km s⁻¹; Geller et al. 2008). Moving at > 10σ_v, these stars are likely on unbound trajectories exceeding the local escape velocity.

7.2 Orbital history reconstruction

Specifically, our two candidates are:

• BSS-E1 (v_rel = 8.9 km s⁻¹, cos θ = 0.99) is currently located at a projected distance of ~14 pc from the core. Given its current velocity, the time since ejection is estimated at τ_ej ≈ 1.5 Myr;

• BSS-E2 (v_rel = 5.9 km s⁻¹, cos θ = 0.98) is also a confirmed non-single star (NSS) candidate, further supporting a dynamical origin involving multi-body interactions.

7.3 The astrometric singularity of BSS-E2: Disruption or a hidden companion

Photometrically, BSS-E2 is firmly located on the red sequence, a demographic conventionally associated with a mass-transfer (MT) formation history within a binary system. However, its Gaia DR3 astrometric parameters present a compelling constraint: with a renormalised unit weight error (RUWE) of 1.00 and a NSS flag of 0, BSS-E2 exhibits no detectable astrometric excess noise, indicating it is currently observed as a single star.

The paradoxical existence of an isolated, high-velocity MT-product BSS points toward two distinct dynamical scenarios:

Scenario I: dynamical ionization. If BSS-E2 was ejected via a violent binary-binary or binary-single encounter, the interaction kinetic energy may have exceeded the internal binding energy of the progenitor binary. In such highly energetic scattering events, the binary system is disrupted (ionized), and the BSS is ejected as a solitary runaway. The observed relative velocity of v_rel = 8.9 km s⁻¹ is significantly higher than the core escape velocity of NGC 188 (v_esc ≈ 3–4 km s⁻¹), confirming that the encounter was sufficiently energetic to severely alter or completely break the orbital configuration of the progenitor.

Scenario II: a hidden white dwarf companion. Alternatively, BSS-E2 may still reside in a post-MT binary system where the companion is a low-mass white dwarf (WD). Since the WD contributes negligibly to the total system luminosity, the astrometric displacement (wobble) depends entirely on the reflex motion of the BSS. At the distance of NGC 188 (d ≈ 1.9 kpc), a BSS (M₁ ≈ 1.5 M_⊙) paired with a typical WD (M₂ ≈ 0.5 M_⊙) in a moderately wide orbit (a ~ 1 AU) induces a maximum astrometric semi-amplitude of merely α ≈ 0.13 mas. Astrometric signature is comparable to or smaller than the single-epoch precision of Gaia DR3. This wobble is perilously close to the single-epoch along-scan precision limit of Gaia DR3, explaining why such a system would seamlessly mimic a single star (RUWE ≈ 1.0).

Regardless of whether BSS-E2 was forcefully ionized or simply harbors an astrometrically invisible companion, its extreme kinematics unambiguously requires a strong dynamical scattering event. This confirms that the MT channel and the dynamical collision channel are not mutually exclusive; rather, dynamically active cluster cores can efficiently scatter MT-formed BSSs into the galactic halo.

Fig. 8 Escape velocity profile of NGC 188 as a function of distance from the cluster center. The dashed red and blue lines mark the current positions of BSS-E1 and BSS-E2, respectively. Both stars possess velocities exceeding the local escape speed, confirming their unbound status.

7.4 Escape velocity and the unbound status

To determine whether the observed runaway BSSs are truly escaping the cluster or merely on high-eccentricity halo orbits, we calculate the local escape velocity profile of NGC 188. Assuming a Plummer potential (Plummer 1911) for the cluster, the escape velocity v_esc at a distance r from the center is given b $$v_{\mathrm{e}\mathrm{s}\mathrm{c}}(r)=\sqrt{\frac{2G{M}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{\sqrt{r^2+a^2}}},$$(8)

where M_tot ≈ 2180 M_⊙ is the total mass derived in Section 3, and a is the Plummer scale parameter, related to the half-mass radius by a = r_h/1.305 ≈ 4.3 pc.

At the cluster core (r → 0), the escape velocity is v_esc,0 ≈ 2.3 km s s⁻¹. At the current projected positions of BSS-E1 and BSS-E2 (r ≈ 12.5 and 21.8 pc, respectively), the local escape velocities drop to ~1.3 km s⁻¹ and ~0.9 km s⁻¹. The observed relative velocities (v_rel = 5.9 km s⁻¹ and 8.9 km s⁻¹) exceed these thresholds by factors of 4 to 10 (see Figure 8).

This significant velocity excess confirms that both BSSs are kinematically unbound from NGC 188. Such high-velocity ejections are a hallmark of strong gravitational scattering involving hard binaries in the core. The fact that these stars are still located within the extended halo despite being unbound is consistent with our orbit back-integration (see next Section): at their current speeds, they would have traveled from the core to their present locations in only ~1.5 Myr, a negligible fraction of the cluster’s 7 Gyr lifespan. These stars thus represent the first direct observational evidence of NGC 188 actively populating the Galactic field with rejuvenated stars.

To further verify the ejection scenario, we reconstructed the past trajectories of the two high-velocity BSS candidates by integrating their orbits backward in time within the cluster potential. We modeled the gravitational field of NGC 188 using a static Plummer potential (Equation (4)), where M_cl = 2180 M_⊙ is the total cluster mass and a = r_h/1.3 ≈ 4.3 pc is the scale radius corresponding to a half-mass radius of r_h = 5.59 pc. We neglect the external Galactic tidal field for this short-term integration (<10 Myr), as the internal dynamics dominate in the cluster core region.

Using the observed positions and relative proper motions converted to physical velocities (assuming a distance of 1.9 kpc), we integrate the equations of motion backwards for 5 Myr using a standard symplectic integrator. The resulting trajectories are shown in Figure 9.

The simulation reveals that:

1. Both candidates trace back directly to the inner region of the cluster (r < r_h) within the last 1.5–2.5 Myr;

2. Specifically, Candidate 1 (ID ...92832) crosses the cluster core at t ≈ −1.5 Myr with a peri-center distance r_peri < 0.5 pc;

3. Candidate 2 (ID ...73024) originates from a similarly dense region at t ≈ −2.1 Myr.

These short timescales (τ_fligh ≪ τ_relax) are consistent with a recent, violent ejection event. The fact that their backward-integrated paths intersect the high-density core – where stellar encounter rates are highest – strongly supports the hypothesis of dynamical formation via binary-binary or binary-single interactions.

Fig. 9 Reconstructed orbital histories of the two high-velocity BSS candidates in the frame of NGC 188. The dashed lines show the backward-integrated trajectories over the last 5 Myr in a Plummer potential (M = 2180 M⊙). The shaded gray circle represents the cluster’s half-mass radius (rh ≈ 5.59 pc). Both stars originate from the deep core of the cluster (r < 1 pc), consistent with a dynamical ejection origin.

7.5 Statistical estimation of the BSS ejection rate

The identification of two high-velocity BSSs with an average ejection age of τ_ej ≈ 1.5 Myr within a single observational snapshot suggests a remarkably active dynamical state for an ancient cluster like NGC 188. Assuming the detected runaway population is representative of the cluster’s current ejection efficiency, we can estimate the time-averaged ejection flux (Φ_BSS) as $${\mathrm{\Phi }}_{\mathrm{B}\mathrm{S}\mathrm{S}}\approx \frac{N_{\mathrm{e}\mathrm{j}}}{{\tau }_{\mathrm{e}\mathrm{j}}}\cdot \frac{1}{C}\cdot f_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{o}},$$(9)

where N_ej = 2 is the observed count, C represents the survey completeness, and f_halo is the geometric correction factor for the unexplored halo regions.

Based on our membership catalog, we derive a mean ejection rate of ~1.33 stars Myr⁻¹ (or ~1330 BSSs Gyr⁻¹). Given the small number statistics, we apply Poisson interval estimation to determine the 95% confidence level. Following the prescription for low-count regimes (Gehrels 1986), the true rate likely falls within the range of [0.16, 4.81] stars Myr⁻¹.

As illustrated in Figure 10, assuming this ejection rate has remained stochastically constant over the cluster’s ~7 Gyr lifespan, NGC 188 would have contributed approximately ${9333}_{-8200}^{+24300}$ stars to the Galactic field. This cumulative loss significantly exceeds the current surviving population (N_mem = 1652), implying that NGC 188 was originally much more massive.

This “rejuvenation engine” effect suggests that old OCs are not merely passive relics; they are active sources of chemically anomalous, high-velocity stars for the Galactic disk. A significant fraction of the so-called “field blue stragglers” or young-looking A-type stars in the solar neighborhood may in fact be “dynamical refugees” ejected from ancient clusters during violent core-scattering events (Perets 2009).

Fig. 10 Statistical prediction of stellar contribution from NGC 188 to the Galactic field. The solid blue line represents the estimated cumulative number of ejected stars over the cluster’s 7 Gyr lifespan, based on a mean ejection rate of ~1.33 stars Myr−1 derived from the two observed runaway BSSs and their ~1.5 Myr ejection age. The shaded light blue region denotes the 95% confidence interval calculated using Poisson statistics for small-sample regimes (Gehrels 1986). For comparison, the horizontal dashed gray line indicates the current census of 1652 confirmed cluster members. While the upper bound is statistically high due to the small Nej, the median estimate of ~9333 robustly illustrates the magnitude of mass loss. This underscores the significant role of ancient OCs in populating the Galactic field with high-velocity, dynamically processed stellar objects.

7.6 Spectroscopic cross-matching and RV limitations

To further characterize the physical properties and 3D kinematics of the two runaway BSS candidates, we performed a cross-match between our catalog and several major spectroscopic surveys, including LAMOST DR10 (Cui et al. 2012), SDSS/APOGEE (Majewski et al. 2017), and the SIMBAD database. Unfortunately, no spectroscopic counterparts were found for these two specific stars in the current data releases.

The lack of spectroscopic data prevents a direct measurement of their radial velocities (V_r), meaning our current dynamical analysis relies on the tangential velocities derived from Gaia proper motions. Consequently, while the extreme tangential alignment (cos θ ≥ 0.98) and high transverse speeds strongly point toward a dynamical ejection from the cluster core, we cannot yet provide a definitive 3D velocity vector or confirm their binary status through multi-epoch RV variations. High-resolution spectroscopic follow-up for these two targets is highly desirable, as it would not only refine their escape trajectories but also reveal atmospheric parameters (T_eff, log g, and [Fe/H]) necessary to distinguish between various BSS formation channels.

8 Conclusions

In this study, we have taken advantage of the unprecedented precision of Gaia DR3 to conduct a forensic kinematic analysis of the old OC NGC 188. Our primary conclusions are as follows:

1. Expanded census and refined photometry: by applying the UPMASK algorithm, we produced a membership catalog of 1652 stars, significantly extending the spatial coverage into the cluster halo. The application of a differential reddening correction reduced the MS width to 0.0424 mag, allowing for the first time a clear photometric bifurcation of the 45 identified BSSs into distinct formation sequences in this cluster;

2. The smoking gun of dynamical ejection: we report the discovery of two high-velocity BSS candidates (BSS-E1 and BSS-E2) in the cluster outskirts. Their extreme kinematics (v_rel ≫ v_esc) and orbital trajectories pointing unambiguously back to the cluster center serve as a kinematic fingerprint of recent dynamical scattering events in the core;

3. Ejection mechanism and nature: the ejection age of ~1.5 Myr suggests that energetic interactions – likely binary-binary or binary-single scattering – are ongoing in NGC 188. While BSS-E2’s position on the red sequence suggests a mass-transfer origin, its current single-star astrometric status (RUWE = 1.0) hints at a violent history involving either dynamical ionization of its progenitor binary or the presence of an invisible white dwarf companion;

4. Impact on the Galactic field: our statistical estimation yields a sustained BSS ejection rate of ~1.33 stars Myr⁻¹. This implies that over its 7 Gyr lifetime, NGC 188 may have deposited nearly 9000 rejuvenated stars into the Galactic disk. This finding suggests that a significant fraction of “field blue stragglers” are likely dynamical refugees from evaporating ancient clusters.

Looking forward, high-resolution spectroscopic observations are essential to confirm the 3D kinematics and binary properties of these runaway candidates. Such data would provide a definitive test for the dynamical ejection model in old OCs and help clarify the role of core-driven processes in the long-term evolution of the Galactic disk population.

In summary, NGC 188 is far from a “quiescent” system; it remains a persistent source of stellar ejections. This work underscores the importance of wide-field kinematic surveys in understanding the full evolutionary impact of OCs on the Galactic stellar population.

Data availability

Python code and data related to the Bayesian modeling algorithm are available on the GitHub repository: https://github.com/chihuanbin/ngc188.

Acknowledgements

This research utilized numpyro, jax, and corner.py. We thank the referee for valuable suggestions. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). This work is supported by the Scientific Research Fund of the Yunnan Provincial Department of Education (No. 2026J0823). This work was also supported by the On-campus Construction Project of the Artificial Intelligence Laboratory, Yunnan Open University.

References

All Tables

All Figures

	Fig. 1 Four-panel comparison between our membership results (blue, N = 1652) and the Hunt & Reffert (2024) catalog (orange, N = 1171). Left: vector Point Diagram in proper motion space. Second: spatial distribution in RA and Dec. Third: histogram of Gaia DR3 parallaxes. Right: G vs. (GBP − GRP) CMD. Our method identifies a significant population of halo members that were previously overlooked.
In the text

	Fig. 2 Differential reddening correction diagnostics. Left: raw CMD with the fiducial MS ridge line (red). Middle: spatial map of the color residuals δE(BP − RP), indicating a fairly uniform extinction field with subtle gradients. Right: final corrected CMD showing sharpened stellar sequences.
In the text

Fig. 3 Structural and dynamical properties of NGC 188. Left: surface density profile derived from the GAIA DR3 membership catalog. The solid line represents the best-fit King model, yielding a core radius, Rc = 2.00 ± 0.15 pc, and a tidal radius, Rt = 51.2 ± 4.3 pc. The background density level, ρbg, is subtracted to highlight the cluster’s extent into the halo. Right: mass segregation ratio, ΛMSR, as a function of the number of the most massive stars, NMST. The significant deviation from unity (ΛMSR > 1.2 for the top 50 massive stars) confirms a strong mass segregation effect, indicating that the cluster is dynamically well relaxed because of its 7 Gyr age.

In the text

Fig. 4 Probability-weighted cumulative mass profile for NGC 188. The solid dark curve represents the cumulative stellar mass M(< R) derived from our corrected catalog (Mtot = 2180 M⊙). The intersection of the dashed red lines identifies the projected half-mass radius rh = 5.59 pc. The shaded region denotes the volume containing 50% of the cluster’s total dynamical mass, providing the structural basis for the adopted Plummer potential.

In the text

	Fig. 5 Color–magnitude diagram of NGC 188 color-coded by estimated stellar mass. The red asterisks mark the identified runaway BSS candidates (see Section 7).
In the text

Fig. 6 Comprehensive analysis of the BSS population. Panel a: CMD showing the BSS selection box relative to the Padova isochrone and the empirical M67 ZAMS. Panel b: histogram of color offsets from the ZAMS, defining the blue and red sequences. Panel c: multiparameter binary identification matrix. Circle sizes are proportional to square of RUWE, and colors represent the image morphology (ipd_frac_multi_peak).

In the text

Fig. 7 Kinematic map of the BSS population in NGC 188. Arrows represent relative proper motion vectors scaled to 8 times their annual displacement. The crimson stars denote the two identified ejected candidates, characterized by high velocities and high alignment scores (cos θ). The core and halo limits are indicated by dashed and dotted circles, respectively. Black crosses denote excluded interlopers.

In the text

	Fig. 8 Escape velocity profile of NGC 188 as a function of distance from the cluster center. The dashed red and blue lines mark the current positions of BSS-E1 and BSS-E2, respectively. Both stars possess velocities exceeding the local escape speed, confirming their unbound status.
In the text

Fig. 9 Reconstructed orbital histories of the two high-velocity BSS candidates in the frame of NGC 188. The dashed lines show the backward-integrated trajectories over the last 5 Myr in a Plummer potential (M = 2180 M⊙). The shaded gray circle represents the cluster’s half-mass radius (rh ≈ 5.59 pc). Both stars originate from the deep core of the cluster (r < 1 pc), consistent with a dynamical ejection origin.

In the text

Fig. 10 Statistical prediction of stellar contribution from NGC 188 to the Galactic field. The solid blue line represents the estimated cumulative number of ejected stars over the cluster’s 7 Gyr lifespan, based on a mean ejection rate of ~1.33 stars Myr−1 derived from the two observed runaway BSSs and their ~1.5 Myr ejection age. The shaded light blue region denotes the 95% confidence interval calculated using Poisson statistics for small-sample regimes (Gehrels 1986). For comparison, the horizontal dashed gray line indicates the current census of 1652 confirmed cluster members. While the upper bound is statistically high due to the small Nej, the median estimate of ~9333 robustly illustrates the magnitude of mass loss. This underscores the significant role of ancient OCs in populating the Galactic field with high-velocity, dynamically processed stellar objects.

In the text