© ESO 2020

1. Introduction

The build-up of the stellar halo of the Milky Way (MW) has been a matter of extensive research for more than five decades. Existing formation scenarios distinguish between “in situ” star formation – where stars form within the host galaxy – and “ex situ” channels, where stars are born in satellite systems and only later accrete onto the massive galaxy. The proposed ratio of these two scenarios spans from purely in situ formation (e.g., Eggen et al. 1962) to mixtures of both channels, with the degree of ex situ contributions varying with Galactocentric distance (e.g., Searle & Zinn 1978; De Lucia & Helmi 2008; Zolotov et al. 2009; Pillepich et al. 2015; Cooper et al. 2015).

There is a wealth of observational evidence of ongoing accretion in the form of persistent stellar streams, which support the ex situ scenarios. These streams range from events as massive as the accretion of the Sagittarius (Sgr) dwarf spheroidal galaxy (Ibata et al. 1994; Belokurov et al. 2006) to observations of more elusive stellar streams and envelopes attributed to globular clusters (GCs) that are in the process of being tidally disrupted (e.g., Odenkirchen et al. 2001; Grillmair & Dionatos 2006; Jordi & Grebel 2010; Kuzma et al. 2018).

The advent of the Gaia space mission (Gaia Collaboration 2016) has revolutionized the field of Galactic archeology. Its latest data release 2 (DR2, Gaia Collaboration 2018a) provides positions, parallaxes, and proper motions for approximately one billion stars (Lindegren et al. 2018). Moreover, for about seven million targets it offers the full 3D space-motion vector by adding radial velocities (Cropper et al. 2018). Gaia DR2 led to the discovery of numerous kinematical substructures comprised of the remnants of massive merger events, such as the incorporation of the Gaia-Enceladus (Helmi et al. 2018; Belokurov et al. 2018) and Sequoia galaxies (Myeong et al. 2019) into the inner Milky Way halo and thick disk (see, e.g., Forbes 2020, for the discussion of potential further merger events), but also many new indications of past and ongoing events of tidal disruption of GCs (Malhan et al. 2018; Ibata et al. 2019; Borsato et al. 2020).

Whether or not a sizable fraction of the halo originated from such disrupted clusters is a long-standing question that remains to be elucidated. One promising approach lies in chemical tagging of halo field stars that show unique chemical signatures in their light-element abundances that are otherwise only found in GCs. This abundance anomaly is attributed to high-temperature proton burning in the CNO-cycle and its Ne-Na chain. These give rise to the archetypal C-N and O-Na abundance anticorrelations found in nearly every GC studied to date (e.g., Cohen 1978; Carretta et al. 2009). The less-frequently observed spreads in Mg, Al, Si, and potentially even K and Zn (e.g., Gratton et al. 2001, 2019; Hanke et al. 2017; Pancino et al. 2017) are attributed to burning at even higher temperatures (≳70 × 10⁶ K, Arnould et al. 1999). Here, we follow the nomenclature of Bastian & Lardo (2018) and distinguish a second, chemically “enriched” GC population (2P), which manifests in enhanced N and Na abundances at the expense of depleted C and O. Models suggest that this second generation was enriched by the aforementioned nucleosynthesis in massive stars of the first, chemically “normal” (i.e., C- and O-enhanced while being Na- and N-depleted) population (1P). While the latter is chemically indistinguishable from the general field star population, discoveries of chemical fingerprints in the halo that are reminiscent of 2P stars are in strong favor of a GC origin.

In turn, a number of theoretical studies investigated aspects of GC formation and evolution by focusing on the amount of mass loss that eventually leads to a direct inference of the fraction of GC stars that contributed to the present-day observed stellar halo (D’Ercole et al. 2008; Bastian et al. 2013; Baumgardt & Sollima 2017; Reina-Campos et al. 2020). Such models ultimately have to be informed by observations of stars from both 1P and 2P, the latter of which can be unambiguously tied to GCs. In this respect, Martell & Grebel (2010) and Martell et al. (2011) used low-resolution spectra of the Sloan Extension for Galactic Understanding and Exploration (SEGUE-1 and SEGUE-2, Yanny et al. 2009; Eisenstein et al. 2011) as part of the Sloan Digital Sky Survey (SDSS, York et al. 2000) to identify 2P candidates in the halo from a combination of CN and CH bandstrengths in red giant branch (RGB) stars. From their finding that 2.5% of their halo sample is CN-strong, these latter authors concluded that between 17 and 50% of the halo may have originated from GCs. Similar arguments hold for an identification of GC-like stars, enabled by the infrared APOGEE survey (Majewski et al. 2016), through their Mg and Al patterns as a product of hotter proton-burning cycles. Respective studies have been successfully carried out for the Galactic halo (Martell et al. 2016; Fernández-Trincado et al. 2019a,b), and others also detected N-enriched stars in the Milky Way bulge (Schiavon et al. 2017), hinting at a similarly (in-)effective formation channel of this old Galactic component.

Using the recent DR14 of SDSS-IV (Abolfathi et al. 2018) we doubled the number of known CN-strong stars by Martell & Grebel (2010) and Martell et al. (2011) to 118 (Koch et al. 2019, hereafter Paper I). From these, we estimated a fraction of 2.6 ± 0.2% 2P stars among all considered halo field giants, which led to a halo contribution from disrupted GCs of 11 ± 1%.

In the present study, we aim to explicitly tie halo field stars – both of enriched and unenriched nature – to a potential GC origin and to observationally test the fraction of bona fide 2P stars among the entire population of GC escapees. To this end, we combined radial velocities (v_r) and metallicities ([Fe/H]) from SDSS/SEGUE with Gaia DR2 astrometry to expand the previous chemical data used in Paper I into a seven-dimensional, chemodynamical information space (i.e., a three-component space vector, a three-component motion vector, and metallicty, [Fe/H]). This was achieved using three methods: First, we searched the immediate vicinity around GCs for extratidal stars that share the space motion and [Fe/H] of the clusters. Next, we used the conservation of action-angle coordinates to identify potential former host clusters for the sample of CN-strong stars from Paper I among the GCs that still exist in the MW. The third method is a modification of the first, and aims at finding stars that share the same portion of phase space as the CN-strong stars, which may indicate a common GC progenitor.

This paper is organized as follows: in Sect. 2, we introduce the data sets employed throughout this work and discuss the involved statistical uncertainties and systematic errors. Next, in Sect. 3, the three employed tagging techniques are outlined, which is followed by the associated findings and a discussion thereof in Sect. 4. Finally, we summarize our main conclusions and provide an outlook for further studies in Sect. 5.

2. Data

For the present analysis, we employed the latest data release of SDSS (DR14, Abolfathi et al. 2018), which contains optical stellar spectra at low resolution (R ∼ 2000) that were obtained throughout the two phases of the SEGUE (Yanny et al. 2009; Eisenstein et al. 2011) survey and during the Extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016). Of particular interest are the large samples of stellar metallicities and heliocentric radial velocities deduced using the SEGUE stellar parameter pipeline (SSPP, Allende Prieto et al. 2008; Lee et al. 2008a,b, 2011; Smolinski et al. 2011).

A large fraction of this study builds upon the compilation of CN-strong giants that were identified in Paper I. There, SDSS spectra of 4649 bona fide halo red giants in the metallicity regime −1.8 dex <  [Fe/H] <  −1.3 dex were examined for a combination of spectral indices involving carbonaceous molecules that are indicative of an unusually strong CN enrichment that is atypical for the general halo population.

In order to make use of the full 6D phase-space information, we complemented SDSS data with proper motions (μ_α cos δ and μ_δ) and parallaxes (ϖ) from Gaia DR2 (Gaia Collaboration 2018a; Lindegren et al. 2018). In addition, the mean phase-space vectors of clusters were retrieved from the published collection by Vasiliev (2019), which itself is based on the compilation of coordinates and distances in the Harris catalog of GC parameters (Harris 1996, 2010 edition) and on line-of-sight velocities from Baumgardt et al. (2019). Other properties used throughout this paper are cluster core radii (r_c), King-model central concentrations (c), and half-light radii (r_h), which were taken from the Harris catalog.

It is important to bear in mind that incorporating both realistic precisions and accuracies is essential for obtaining realistic probabilities in our Bayesian framework (Sect. 3). This latter is used to compare quantities from different sources and is therefore affected by potential systematic discrepancies. Below, we discuss the two main sources of error that were singled out throughout our analysis, namely inaccurate metallicities with additionally underestimated precisions, and underestimated errors for Gaia proper motions.

2.1. SSPP parameters

In light of the low-resolution (R  ∼  2000) nature of the SEGUE/eBOSS surveys and possible systematic trends among the parameter scales involved in the SSPP, caution must be taken when using the quoted internal uncertainties on the adopted parameters. A realistic assessment of the error budget should incorporate both systematic error contributions such as wavelength calibration errors or calibration inaccuracies of the SSPP parameter scales, as well as random errors caused by the finite signal-to-noise ratios (S/Ns) of the underlying spectra. Substantial metallicity- and temperature-dependent departures of the SSPP stellar parameter scales with respect to high-resolution reference samples have already been pointed out by Smolinski et al. (2011; see, in particular, their Fig. A2). The authors further compared SSPP results for [Fe/H] and v_r for several GCs with literature values and claim overall “good internal and external consistency” despite the evidence they provide for a wealth of residual substructure. However, an independent analysis of SEGUE spectra using the flux-ratio-based parameterization method ATHOS (Hanke et al. 2018) shows very similar residual trends with respect to SSPP. This observation contradicts the favored explanation by Smolinski et al. (2011), who speculate that deviations may originate from inhomogeneities among the high-resolution studies that were used as reference.

While performing an in-depth evaluation of different error sources, we compiled a sample of stars with a high probability of being associated with the GC M 13. This cluster was chosen because it is well-studied and, due to its proximity, has been targeted with many SDSS fibers. Only SDSS targets that fall within a projected separation of one tidal radius from the cluster center were considered. Furthermore, we employed Gaia parallaxes and v_r from SDSS as a means to reject fore- and background stars by selecting only those with insignificant (< 5σ) deviations from the mean values of the cluster. We note that at this point neither information about the chemical composition from SDSS nor Gaia proper motions entered the selection procedure. Nevertheless, as can be seen in Fig. 1, the resulting color–magnitude diagram (CMD) based on Gaia DR2 G and G_BP − G_RP photometry appears clean with a low degree of contamination (likely less than five stars out of 283). On this account, the vast majority of the sample can be assumed to be cluster members. When comparing SSPP results for [Fe/H] with the literature value (−1.53 dex, Harris 1996), systematic discrepancies as a function of the evolutionary state of the individual stars (Fig. 1) become apparent. Specifically, stars on the upper main sequence (MS) deviate by −0.18 dex with a decreasing trend toward the subgiant branch (Δ[Fe/H] ≈ −0.10 dex), while the RGB is consistent with a zero difference. Stars on the horizontal branch (HB) on the other hand differ by about −0.20 dex. Given the fact that M 13 shows no signs of intrinsic metallicity spread (e.g., Johnson & Pilachowski 2012), any such difference can be attributed to SSPP inaccuracies.

Fig. 1. Precision and accuracy validation of [Fe/H]SSPP for M 13 member stars. Lower left panel: proper-motion-cleaned CMD in the GaiaG, GBP, and GRP passbands. For reference to the lower right panel, GBP − GRP is additionally indicated by colored symbols. Right panels: marginalized distribution of deviations (top) and behavior with evolutionary state (bottom) of SSPP metallicities with respect to the literature value: [Fe/H]M 13 = −1.53 dex. Red points and error bars denote the means and standard deviations in bins of 1 mag width in G, whereas underlying gray bars represent the median of the SSPP uncertainties on [Fe/H].

Again taking advantage of the fact that there are effectively no detectable metallicity spreads in M 13, we can estimate the internal SSPP precisions for [Fe/H] at different evolutionary states from the observed scatter. Typically, this scatter is larger by a factor of three compared to the provided SSPP uncertainties. The latter comparison is illustrated in the right panel of Fig. 1.

Given the caveats discussed above, we quadratically added a global systematic error of 0.15 dex to individual errors on [Fe/H] measurements¹ throughout this work. This drastically increases the majority of the quoted (internal) SSPP uncertainties. We further caution that the total error budget on v_r cannot be less than 5 km s⁻¹ (see the comparison of SSPP radial velocities to high-resolution results by Smolinski et al. 2011), which we accordingly adopted as another systematic error that is added in quadrature.

2.2. Gaia astrometric solution

For stars with an available five-parameter astrometric solution

$$\begin{array}{c}\hfill \mathit{x}={(\alpha ,\delta ,\varpi ,{\mu }_{\alpha }cos\delta ,{\mu }_{\delta })}^{\mathrm{T}},\end{array}$$(1)

the Gaia DR2 data tables enable the computation of the full covariance matrix, Cov(x), for the solution of each individual star. Incorporating its off-diagonal entries is crucial for this work, in that the measurements of ϖ, μ_α cos δ, and μ_δ may be correlated to varying degrees and thus must not be considered independent. Furthermore, we follow the technical note GAIA-C3-TN-LU-LL-124-01² and scale the covariance matrices by the squared re-normalized unit weight error (RUWE),

$$\begin{array}{c}\hfill {\mathrm{RUWE}}^2=\frac{{\chi }^2}{(N-5)·u_0^2(G,G_{\mathrm{BP}}-G_{\mathrm{RP}})},\end{array}$$(2)

where χ² is the chi-square value of the astrometric fit to all N Gaia measurements in the direction along the scan that are considered “good”. The factor u₀ is an empirically calibrated quantity that can be extracted from dedicated lookup tables as a function of G and G_BP − G_RP³. We emphasize that all employed covariance matrices involving Gaia data were scaled by RUWE². In this work we moreover corrected for the quasar-based ϖ zero point of −0.029 mas (Lindegren et al. 2018) and assumed an additional global systematic error for proper motions of 0.035 mas yr⁻¹ (Gaia Collaboration 2018b).

In elaborating on the formal DR2 proper-motion errors we again employ the sample of M 13 stars introduced earlier in this chapter. Their spatial distribution is presented in Fig. 2 next to the corresponding CMD. We further show the distribution of the selection in absolute and relative (i.e., scaled by standard errors) proper-motion differences with respect to the mean values for M 13. While bright stars are highly clustered in an absolute sense, from Fig. 2 it becomes apparent that there is an inversion of the relative distribution widths of bright and faint targets when going from absolute to relative proper motion differences. It is noteworthy that the RUWE remains below 1.5 for all but nine targets in the selection, four of which are fainter than G = 17 mag. Hence, for the vast majority of our bright stars, inseparable blends at the edge of the detection limit can be excluded with high confidence. We therefore attribute our observation to an additional, hitherto unexplained systematic error component of the proper motions in excess of the already applied 0.035 mas yr⁻¹.

Fig. 2. Upper left: same CMD as in Fig. 1, but with color coding indicating the brightness (see lower right color bar); the same color code is used in all other panels. Upper right: spatial distribution of the M 13 validation sample. The cluster core, half-light radius (1.69 arcmin), and tidal radius (21 arcmin) are represented by a black star, red circle, and black-dashed line, respectively. We emphasize that the tidal radius together with ϖ and vr were the only selection criteria applied to investigate proper motion systematic errors (main text). Middle panels: absolute (left) and relative (right) deviation in proper motion of the sample from the cluster mean value. Lower panel: squared Mahalanobis distance, M2, from proper motions only versus projected distance from the cluster center. In analogy to the 1σ, 2σ, and 3σ significances of a normal distribution, dash-dotted, dashed, and dotted lines represent the corresponding p-values (see legend). The red line denotes the scaling relation introduced in Eq. (4).

Although we neglected covariances for the middle panels of Fig. 2 for illustrative purposes, owing to the discussed correlations among the components of y = (μ_α cos δ,  μ_δ)^T, it is not straightforward to estimate the statistical significance of the deviation from the cluster mean, y_GC. Hence, for each star, we computed the squared Mahalanobis distance

$$\begin{array}{c}\hfill M^2={(\mathit{y}-{\mathit{y}}_{\mathrm{GC}})}^{\mathrm{T}}\mathrm{Cov}{(\mathit{y}-{\mathit{y}}_{\mathrm{GC}})}^{-1}(\mathit{y}-{\mathit{y}}_{\mathrm{GC}}),\end{array}$$(3)

which respects the combined covariance. Under the assumption of normally distributed errors, the latter quantity is chi-square distributed. Therefore, the p-value – that is, the probability of finding a value of M² or more extreme under the null hypothesis that the star is not kinematically distinct from the cluster – can be computed from a chi-square distribution with one degree of freedom. In Fig. 2, M² is plotted against the projected angular distance to the cluster center, r, normalized by the half-light radius, r_h. We identify two distinct groups of highly significant outliers that would be rejected in Sect. 4.1 based on their proper motions alone. First, starting from r/r_h ≈ 6, there is a trend of increasingly significant deviation with decreasing r. Secondly, irrespective of the separation from the cluster core, bright stars (G ≲ 16 mag) tend to deviate more significantly than fainter ones. Based on the clean CMD and the v_r-based association, we exclude the possibility that the majority of the strong proper motion outliers do not in fact belong to the cluster and conclude that the origin is genuinely to be found in underestimated Gaia errors. We compensate for the latter by introducing a distance-dependent factor

$$\begin{array}{c}\hfill \alpha (r/r_{\mathrm{h}})=\{\begin{array}{cc}{10}^{-0.3(r/r_{\mathrm{h}}-6)},\hfill & \mathrm{if}r/r_{\mathrm{h}}\le 6\hfill \\ 1,\hfill & \mathrm{otherwise},\hfill \end{array}\end{array}$$(4)

to be applied to the covariances. Though optimized for the M 13 stars, we note that we found the relation to hold true for several other GCs from our sample within an acceptable amount of variance.

3. Methods

In the following sections, we introduce the three approaches that were used to identify potential former GC members. The first method is tailored to test for stellar associations in the immediate cluster surroundings and does not require all of the three space motions and positions to be measured with the same precision. The second approach ties field CN-strong stars to clusters and bears the advantage of not being limited to a comparatively small patch of the sky, though at the expense of obtaining results that are strongly impacted by the lowest-precision entry in the phase-space vector. The third and last method is an adaptation of the first method, where we search for chemically normal stellar populations that share the same metallicity and kinematic properties as CN-strong stars in the halo field, possibly indicating a common origin.

In all three approaches we employed variations of the Bayesian approach described in Anderson et al. (2013). In brief, the association probability (posterior) of a star being attributable to a reference object (be it a CN-strong star or a GC) is given by

$$\begin{array}{c}\hfill P(A|B)=\frac{P(B|A)\times P(A)}{P(B)}\propto P(B|A)\times P(A).\end{array}$$(5)

Here, P(B|A) is the conditional probability (likelihood) of encountering the star at its observed position in information space assuming it is indeed associated with (or used to be an intratidal member of) the reference. Further, P(B) is the probability of observing the data and P(A) is the initial degree of belief in association (prior).

3.1. Method I: Stars in the neighborhood of clusters

The potential GC origin of individual stars in the cluster vicinity (out to a few degrees separation) was evaluated on a cluster-by-cluster basis. Stars that escaped the GC potential and show a significant spatial separation from their previous host – while conserving their metallicity – no longer necessarily share the same space-motion vector as the cluster. We instead presume that their currently observed motion should show a stream-like behavior and thus be similar (though not exactly identical; see Sects. 4.2 and 4.4) to the closest point along the hypothesized stream.

To implement this behavior, we extrapolated the orbits of all GCs in the list by Vasiliev (2019) to 2 Gyr in the past and future by performing point-particle integrations employing the python library galpy (Bovy 2015) and its standard Galactic potential MWPotential2014. As a means to track the behavior of stars lost from the clusters, we further used galpy for the modeling of dynamically cold (σ_v = 1 km s⁻¹) leading and trailing tidal streams. We introduced tracer particles along the streams with parameter covariances that account for the simulated distribution functions of the streams. From these, the likelihood terms in Eq. (5) were calculated individually for each star based on the difference Δz between the five-component vector z = ([Fe/H], v_r, ϖ, μ_α cos δ, μ_δ)^T and the corresponding closest stream anchor point, z_s. Again following Anderson et al. (2013), in analogy to Eq. (3), we computed the squared Mahalanobis distance M². Here we used the Gaia covariance entries for ϖ, μ_α cos δ, and μ_δ and assumed no correlation between SDSS and Gaia quantities. Furthermore, covariances between cluster mean proper motions were taken from Vasiliev (2019). Complete independence of the cluster quantities from the stellar quantities was presumed, such that the cross-covariances between z and z_s are zero. With respect to the second assumption, while our potential extratidal candidates were not used by Vasiliev (2019) to constrain mean cluster parallaxes and proper motions, this does not necessarily hold for our supposedly bound cluster members. Vasiliev (2019) on the other hand commonly employed several orders of magnitude more stars from the Gaia tables than there are counterparts in the SDSS catalog. As a consequence, the overlap and therefore cross-covariances are minor. Interdependencies between ϖ, [Fe/H], v_r, and their respective mean cluster values can be excluded, as the cluster parallaxes (by means of inverse distances) and [Fe/H] are based on the Harris catalog of GC parameters (Harris 1996), and v_r stems from the collection of ground-based measurements compiled by Baumgardt et al. (2019). None of the former are in any way connected to SDSS or Gaia.

Finally, the likelihood can be expressed as

$$\begin{array}{c}\hfill P(B|A)=1-p(M^2),\end{array}$$(6)

where p(M²) is the p-value of a χ² distribution with five degrees of freedom. Anderson et al. (2013) emphasize that large errors – translating into low-significance values – cannot rigidly exclude a large portion of their corresponding phase-space dimensions and thus do not limit the high-likelihood regime to a confined range. A prime example of this behavior is ϖ, which – in light of typically large cluster distances – hardly exceeds the 2σ significance level for most of the stars that are deemed cluster members below. Nonetheless, ϖ is a powerful means to reject the large number of foreground stars, exposing significant parallaxes that are inconsistent with cluster association.

In contrast to Anderson et al. (2013), who demanded their targets be gravitationally bound to the clusters, we cannot use this criterion for extratidal stars and hence assume a loose prior of the form

$$\begin{array}{c}\hfill P(A)=\exp (-\frac{1}{2}{\left(\frac{d_{\mathrm{s}}}{7r_{\mathrm{t}}}\right)}^2),\end{array}$$(7)

with d_s being the projected distance of a star to the closest tracer stream particle and r_t being the tidal radius of the cluster. This prior formulates the initial belief that former members at arbitrarily large angular separation will not share the same z.

Figure 3 illustrates a graphical representation of the multivariate association procedure for the exemplary case of the cluster M 13. Candidates are reported if their associated P(A|B) exceeds the threshold of 0.05. Both intra- and extratidal stars are treated in the same way.

Fig. 3. Five-parameter chemodynamical association criteria used to constrain the association with the GC M 13. The rejected back- and foreground population of stars with SSPP parameters and Gaia DR2 motions is indicated with gray density contours. Brown dots indicate the simulated stream (see main text). Intratidal and extratidal stars are depicted with blue and colored circles (see also Fig. 7), respectively. Circle sizes indicate the association probability and colors resemble the projected distance to the cluster center (see legend and color bar). The marginalized distributions of each coordinate on the abscissa are shown on top of the two-dimensional representations. The adopted spatial prior as well as the stars that are rejected due to its usage (gray circles with no colored counterpart) are illustrated in the upper right-hand panel.

3.2. Method II: Integrals of motion

It is possible that GC escapees that did not recently become unbound from their host cluster may no longer be found in the immediate cluster vicinity; such GCs would not be recovered by the approach in the previous section. Thus, for the manageably small sample of CN-strong stars used here, we resort to the fact that for most orbits in axisymmetric potentials there exists a set of three conservative integrals of motion (e.g., Henon & Heiles 1964); in other words, despite being spatially separated, escapees almost completely retain the integrals of motion⁴ of their host (e.g., Savino & Posti 2019). We used the galpy implementation of the Stäckel approximation by Binney (2012) in order to integrate the axisymmetric MWPotential2014 for the radial and vertical actions, J_r and J_z, as well as the azimuthal component of the angular momentum, L_z. With respect to the previous approach, such integrations bear the main disadvantage of relying on the full six-component phase-space vector to initialize an orbit, such that uncertainties are strongly driven by the least constrained observational quantity. While the latter restriction is not a major concern for our sample of GCs, the lack of highly significant measurements for the heliocentric distance, D, to the CN-strong stars ultimately drives the confidence for rejecting or accepting a potential cluster association. Gaia parallax significances ϖ/σ_ϖ <  1 have been attributed to 61% (68 stars) of our sample whereas only three measurements exceed the 4σ level. It is evident that any distance inferred from ϖ alone (e.g., Bailer-Jones et al. 2018) will result in critical errors in the actions.

In order to obtain more precise estimates for D, we refined the spectrophotometric formalism from Paper I by adopting the Bayesian inference method for stellar distances by Burnett & Binney (2010) using photometry and stellar parameters from SDSS and ϖ as an additional constraint (cf., e.g., Savino & Posti 2019). Therefore, a likelihood was computed from the residuals between the observed quantities (color, apparent magnitude, T_eff, [Fe/H], log g, and ϖ) and their theoretical counterparts from a grid of PARSEC (Marigo et al. 2017) isochrones. As prior, we used the three-component Galactic model adopted by Burnett & Binney (2010) for the thin and thick disks and the stellar halo, the latter being constructed from the parameters stellar age, [Fe/H], and mass. The magnitude and variance of D were deduced by means of the first and second moments of the star’s posterior probability density function (pdf) marginalized in all other dimensions (cf., e.g., Savino & Posti 2019). Even though the nominal errors on the inferred distances are small and consequently the statistical significances are high (nearly 50% reside above D/σ_D = 10), we caution that there is ample room for various systematic errors that are not captured by the mentioned treatment. One example is the SSPP surface gravity – a quantity that is notoriously hard to obtain with any accuracy from low-resolution spectra – which may strongly favor much larger distances in cases where a star was erroneously classified as a low-gravity giant whilst in fact being a giant of intrinsically higher gravity or even a dwarf star.

Adopting the former distances, we present the resulting spatial distribution of our CN-strong giants in the Galactocentric frame on the right-hand side of Fig. 4, while a Hammer projection in the Galactic frame is shown on the left-hand side. Furthermore, in Fig. 5, we show a Toomre diagram of our CN-strong stars and GCs of the same [Fe/H], from which, qualitatively, a remarkable distribution overlap can be seen. Nonetheless, we envision a more quantitative approach, exploiting the kinematic memory of our targets.

Fig. 4. Spatial distribution of GCs (blue squares) and CN-strong halo giants (colored according to their distance) in Galactic coordinates (left) and the Cartesian Galactocentric frame (right). The dashed line in the left panel denotes the celestial equator. The SDSS footprint area is depicted in white while the Gaia all-sky map is illustrated in the background for orientation. Image credit. Gaia Data Processing and Analysis Consortium (DPAC); A. Moitinho/A. F. Silva/M. Barros/C. Barata, University of Lisbon, Portugal; H. Savietto, Fork Research, Portugal.

Fig. 5. Toomre diagram for our CN-strong giants (colored circles) and the Galactic GC population (blue squares). The latter is only shown in the regime −1.8 <  [Fe/H] <  −1.3 dex, representing the coverage of the analyzed stars. The color coding of the CN-strong stars is the same as in Fig. 4. For better comparability, an additional Gaussian kernel density estimate of the GC sample is provided where increasing color darkness resembles increasing density. To show the clear contrast to the Galactic disk population, we also indicate all targets from our Gaia – SDSS crossmatch that obey −5° < b <  5° and have a spectrophotometric distance as determined by SDSS of less than 1 kpc (gray dots). Moreover, the dashed line encircles the region of absolute velocities of less than 100 km s−1 with respect to the local standard of rest.

To this end, all a priori known statistical error sources on the phase space vectors for both GCs and CN-strong stars were propagated in the action integration by means of a Monte Carlo analysis. Each orbit was initialized 500 times with values randomly drawn from a multivariate normal distribution respecting the covariance matrix, where we assumed cross-covariances between α, δ, D, and v_r to be negligible while maintaining the covariances for Gaia proper motions. As shown by the exemplary comparison in Fig. 6, the obtained distributions in action space are highly coupled and non-normal, that is, the relations between action coordinates show a strong nonlinear behavior. It is apparent that, whilst cluster (M 70) and star (Gaia DR2 3833963854548409344) indicate a perfect match in the marginalized, one-dimensional distributions, in the multidimensional representation there is only a much smaller distribution overlap. If not properly accounted for, this effect would overestimate the association probability in most cases.

Fig. 6. Distribution of Monte Carlo samples for the integrals of motion for M 70 (red) and the star Gaia DR2 3833963854548409344 (blue). The marginalized two-dimensional and one-dimensional kernel density estimates are indicated as contours in the frames and histograms at the top, respectively. For illustrative purposes, the histograms were normalized by their maxima.

To test the possibility of former membership of a CN-strong halo star to a cluster, we required both the actions and [Fe/H] to be consistent. Hence, we compared the cluster vector ${(J_{\mathrm{r}},J_z,L_z,[\mathrm{Fe}/\mathrm{H}])}_{\mathrm{GC}}^{\mathrm{T}}$ to the corresponding vector for the star. Due to the non-normality of the MC samples mentioned above, we cannot simply compute the covariance matrices of the samples and proceed with the Mahalanobis distance employed in Sect. 4.1. Therefore, we approximated the four-dimensional probability density function, ψ(v), of the difference vector, v, using a nonparametric Gaussian kernel density estimate. We note that this is the main difference between our approach and that of Savino & Posti (2019), who did not account for nonlinear correlations between actions. As a means to accept or reject association, we obtained

$$\begin{array}{c}\hfill p={\displaystyle \underset{\psi ({\mathit{v}}^{\prime })<\psi (\mathbf{0})}{\int }\psi ({\mathit{v}}^{\prime })\mathrm{d}{\mathit{v}}^{\prime }}\end{array}$$(8)

through Monte Carlo integration.

3.3. Method III: Chemodynamical matches in the fields surrounding CN-strong stars

Following the predictions made by evolutionary GC models (e.g., D’Ercole et al. 2008; Bastian et al. 2013), a substantial number of first-generation stars are lost early on in the cluster formation phase. This unenriched population is chemically no different from the standard collection of halo stars and thus remains unidentifiable for classical tagging methods (e.g., Paper I). Yet, by exploring chemodynamical similarities between the enriched, bona fide second-generation cluster stars in the Galactic halo and their surrounding field, we can search for links that could indicate a common origin. Therefore, we once again employed the formalism of Sect. 3.1 with CN-strong stars taking over the role of the GCs. The main difference between this latter method and the one described here is that we did not attempt to integrate the orbit of the CN-strong stars because the much larger errors on the astrometry of these individual stars do not allow for the computation of the tightly constrained orbits required to confidently reject field interlopers. As a spatial prior, in analogy to Eq. (7), we adopted a Gaussian centered on the CN-strong star with a standard deviation of seven degrees in projected distance.

4. Results and discussion

4.1. Extratidal escapee candidates around clusters

In this section, we present the discovery of new candidates for extratidal stars around GCs in the halo field based on their position in multi-dimensional information space. We started by cleaning the SSPP catalog for stars lacking a [Fe/H] or v_r detection. The remaining list was subsequently cross-matched with the Gaia DR2 source catalog. This vetting process – due to the SDSS sky coverage being limited to the northern sky (cf., Fig. 4) – naturally rejected the southern GCs and left us with about 3.7 × 10⁵ stars for further analysis. The relevant cluster properties entering the analysis can be found in Table 1.

We consider those stars with P(A|B) > 0.05 as bona fide cluster associates. A further sub-classification into extra- and intratidal candidates is performed based on whether the projected distances of the stars do or do not exceed the tidal radius of the cluster, r_t. We note the major caveat of this distinction to be the potential exclusion of truly extratidal stars that, owing to their line-of-sight position, happen to be projected into the tidal radius of the cluster. Unfortunately, the already mentioned low significance of ϖ for these distant objects cannot yield additional constraints and prohibits their identification.

In the following sections, we comment on our 151 extratidal candidates with likely associations with eight GCs. The discussion is separated by individual host clusters. A summary of the properties of all halo stars with potential GC origin is provided in Table C.1.

4.1.1. M 13 (NGC 6205)

M 13 is an old (12.3 Gyr, VandenBerg et al. 2013), intermediately metal-poor (−1.5 dex, Johnson & Pilachowski 2012) GC known to host extreme light-element abundance variations (e.g., Sneden et al. 2004; Johnson & Pilachowski 2012). It has been shown by Cordero et al. (2017) that the most strongly enriched population of stars possesses a stronger degree of rotation with respect to the other samples, while Savino et al. (2018) were not able to find significant display of spatial segregations. Jordi & Grebel (2010) reported on an extratidal stellar halo in the immediate vicinity of this GC with a slight elongation aligned with its motion that was later extended to even larger distances out to 13.8 r_t by Navin et al. (2016; however, we note the reevaluation of these associations later in this section).

For M 13, we identified a total of 292 candidates that – according to our classification scheme – are candidate cluster associates. Out of those, 260 are still bound to the cluster and 32 are most likely extratidal. A multi-parametric representation of the chemodynamical associations against the rejected background population is presented in Fig. 3. The clustering of the associates around the mean values for M 13 and the simulated stream is evident. Interestingly, M 92 (Sect. 4.1.2) can be identified as a slight background overdensity. Nonetheless, considering the full five-dimensional information space, M 13 can be discerned from M 92 stars with high confidence.

An independent and powerful constraint for the validity of our method can be obtained from the color–magnitude information for our candidates, which is presented in Fig. 7 along with the spatial distribution. Though they are not part of the initial analysis, the distributions in the Gaia and SDSS passbands constitute well-behaved structures, which – with few exceptions and considering observational uncertainties – are consistent with a single ∼12 Gyr PARSEC isochrone that matches the age, [Fe/H], and distance modulus of M 13 (see Table 1). We note that the reddening value for the Gaia isochrone displayed in Fig. 7 had to be manually adjusted to achieve a better fit. When using reddening coefficients provided in Yuan et al. (2013) for SDSS colors and from Casagrande & VandenBerg (2018) for Gaia colors, this reddening is not consistently described by the same visual extinction A(V) as the SDSS colors. A detailed investigation of this behavior is beyond the scope of the present study. Due to this and the effects of varying extinction throughout the substantial sky coverage of more than 10 deg away from the cluster center, we do not attempt to incorporate the photometry in our mathematical treatment, but use it to qualitatively assess membership likelihoods. To guide the eye in the Gaia CMD, we furthermore show all Gaia targets – not just the ones with a counterpart in SDSS – that lie within the tidal radius and that do not deviate by more than 3σ in ϖ, μ_α cos δ, and μ_δ. Additionally, in order to overcome contamination by stars that might be affected by crowding and therefore might be blueshifted from the MS, the MS turnoff (MSTO), and the RGB, we required the Gaia color excess factor⁵ to remain below 1.4.

Fig. 7. Parameter distributions of intra- (blue) and extratidal (black to light yellow) stars associated with M 13 on top of the field star population covered by the SDSS footprint (gray dots). The sizes of the colored circles are directly proportional to the membership probability P(A|B) as indicated in the legend, while the color of the extratidal candidates denotes the projected distance from the cluster core in units of rt (see color bar). Upper panels: CMDs from Gaia (left) and dereddened SDSS (right) bands. Cyan dots indicate proper motion and ϖ-selected Gaia sources without spectroscopic SDSS counterparts (see main text for details). Typical photometric errors are specified in each case by error bars on the left-hand side. For a subset of bright intratidal stars, reliable color information from SDSS is not available, which is why those targets were omitted in the right CMD. Red lines are 12 Gyr PARSEC isochrones (Marigo et al. 2017) matching the parameters from Table 1 with an applied manual shift in color for the Gaia CMD. The CN-strong giants classified in Paper I are labeled “CN” in both panels. Lower panels: distribution on the sky (left) and a zoomed-in view of the immediate cluster vicinity (right). The black-dashed ellipse indicates the tidal radius, rt, and the direction to the Galactic center is marked by a gray arrow. Properties for the integrated cluster orbit and the simulated stream are indicated by the blue curve and brown dots, respectively.

The projected angular distances of 13 out of the 32 extratidal candidates are only slightly larger than the tidal radius (d <  1.2 r_t) and thus may still be considered loosely bound to M 13. Their distribution qualitatively matches the photometrically determined contours by Jordi & Grebel (2010). At d = 1.6 r_t and above, the remaining 19 candidates can safely be denoted extratidal. Looking at their spatial distribution (Fig. 7), it is tempting to claim evidence for an overdensity towards the leading portion of the cluster orbit and/or stream. However, owing to the nonisotropic coverage of the SDSS footprint around the cluster, this observation should be treated with caution.

Apart from one star on the blue HB, the three brightest extratidal associations are affiliated with the RGB and were classified as CN-strong in Paper I. Furthermore, the stars are also attributable to M 13 by the method outlined in Sect. 4.2. Therefore, it appears very likely that the three stars in question are escapees that can be attributed to the second – or enhanced – population of M 13 stars. We emphasize that among all 32 candidates presented here, only four giants were studied in Paper I due to strict parameter limitations of the method therein. In other words, 75% of the stars in the region of parameter overlap are CN-strong. This being the case, it is intriguing to argue that a substantial fraction of all M 13 escapees might be second-generation stars, but due to low-number statistics, this remains rather speculative (see Sect. 4.5 for further discussion).

Navin et al. (2016) performed a pre-Gaia-era search for cluster associations both within and outside of the tidal radius of M 13. The study was based solely on space motions and photometry. The authors obtained initial membership estimates relying on radial velocities from the first data release of the LAMOST spectroscopic survey (Luo et al. 2015) and subsequently refined the selection criteria using cuts in color–magnitude and stellar-parameter space in order to end up with probable cluster giants and to exclude foreground dwarf stars. A further constraint was established by employing UCAC4 (Zacharias et al. 2013) proper motions and demanding deviations from the cluster means of less than 10 mas yr⁻¹. We note that the proper motions used in this latter study are barely significant in any of the cases. Due to large uncertainties in LAMOST DR1, [Fe/H] was not used by the authors to further eliminate candidates.

Here, we cross-matched the 12 candidate extratidal halo stars for M 13 from Navin et al. (2016) with the Gaia DR2 catalog (see Fig. 8). Given the bright nature of the stars in question, corresponding Gaia proper motions are highly significant⁶ and from a visual inspection of the proper motion comparison in Fig. 8 it is already evident that most claimed candidates can be excluded with high confidence. Moreover, our membership formalism (Eqs. (5) and (6)) provides a zero probability and thus rejects all 12 candidates. Here, we assumed the sample mean and median errors quoted by Navin et al. (2016) on the LAMOST quantities of 17.1 km s⁻¹ and 0.86 dex for v_r and [Fe/H], respectively. The latter values are relatively conservative, such that our strong exclusion confidence is almost entirely driven by the high-precision Gaia proper motions, which could not even remotely be equalled by available catalogs at the time of publication of the original study. At statistically significant parallax-based radial distances of 4.1 kpc and 4.7 kpc – as opposed to a cluster distance of 7.1 kpc – two stars can be ascribed to the foreground population and are excluded by this parameter alone.

Fig. 8. Comparison of the chemodynamical (left two panels) and spatial (third panel) distributions and the CMD (right) for M 13 between Navin et al. (2016; cyan stars) and the present study (same color and size scheme as in Fig. 7). In the proper motion diagram, the values from the original study are connected by red lines with their Gaia DR2 counterparts (red stars).

4.1.2. M 92 (NGC 6341)

At [Fe/H] = −2.35 dex (Carretta et al. 2009), the GC M 92 is amongst the most metal-poor Galactic clusters known. Using photometric data, Testa et al. (2000) and Jordi & Grebel (2010) consistently reported on an extratidal halo for M 92. It is noteworthy that our debris model for this cluster is exceptionally scattered because of its rather close pericentric passages (R_peri ≈ 0.17 kpc) that result in a strong tidal field.

Figure 9 indicates the distribution in parameter space of our 110 potential M 92 associates, of which 35 are likely to be extratidal. Due to the low metallicity of the cluster, and its considerable proper motion in α cos δ direction (−4.9 mas yr⁻¹), it is straightforward to chemodynamically discern bona fide cluster associates from the fore- and background field stars and hence effectively diminish the false-positive detection probability. From the CMDs we find that – with the exception of one low-probability case at large angular separation – our candidates are consistent with the cluster isochrone and that they are distributed between the upper MS (G ∼ 19 mag) and the upper RGB (G ∼ 14.5 mag). Two stars fall on top of the blue HB of M 92. While nine candidates are borderline cases in terms of their distance to the cluster core (d <  1.30 r_t), 26 stars, mostly positioned south of the cluster, are clearly unbound at core separations of more than 1.3 tidal radii. Unfortunately, M 92 was placed on the northern edge of the SDSS plate such that any potential extratidal candidates in the immediate surroundings north of the cluster are missed by our analysis. The eastward elongation of the distribution of chemodynamically associated stars in the near-field surrounding the cluster shows a striking similarity to the photometric overdensity reported by Jordi & Grebel (2010).

Fig. 9. Same as Fig. 7, but for M 92.

Since our approach is density insensitive and merely limited by the SDSS selection function, we can expand the search for cluster members out to larger separations. Moreover, omitting the clear contaminant from the CMD inspection, we associate 12 targets at 12.5 <  d/r_t <  75.2 with M 92, thus expanding the previously known extension of the tidal debris from this GC.

4.1.3. M 3 (NGC 5272)

Of the 199 associates for M 3, we report on 17 probable candidates in the surrounding field star population that coincide with evolutionary stages between the cluster’s upper MS, base of the RGB, and HB. The brightest extratidal star (G = 15.8 mag) shares the same region as the locus of intratidal HB stars. In fact, both Abbas et al. (2014) and Clementini et al. (2019) list this star as being a pulsating RR Lyrae star of type ab. The former authors provide a distance estimate of 10.99 kpc, which – lacking a proper attached error margin on that value – is hard to compare with the cluster distance of 10.2 kpc, although an association seems feasible. We performed our own analysis of the photometric time series available in the literature and furthermore corrected for the fact that, due to their pulsating nature, RR Lyrae do not straightforwardly reveal their systemic velocity from single-epoch spectra alone. A detailed description can be found in Appendix A. Our inferred value of 10.9 ± 2.6 kpc is in excellent agreement with the cluster distance. Kundu et al. (2019) present another two RR Lyrae stars that they associate with M 3. However, the search radius of these latter authors was restricted to 2/3 <  d/r_t <  3, a circumstance that prohibited them from finding the variable presented here at d = 9.4 r_t.

Among the presumably bound and unbound candidate stars, several reside blueward of the MSTO and subgiant branch (SGB) as can be seen in Fig. 10. These could either be connected to false-positive associations or be genuine blue straggler stars (BSS). We favor the latter interpretation, as we do not see any particular reason for all random associations to be preferentially found on the blue- as opposed to the red side of the isochrone, even though the field star population is much more numerous in the red part. Our explanation is bolstered by the fact that the number of true BSS in M 3 is much larger than in M 13 (Ferraro et al. 2003) for example, where our treatment did not associate any potential blue straggler candidate.

Fig. 10. Same as Fig. 7, but for M 3. The approximate color range of the instability strip from Clementini et al. (2019) is indicated in the upper left panel by vertical dashed lines and the pulsating variable (see main text) is labeled “RRab”. We note that the illustrated color and brightness for the RR Lyrae star are merely mean quantities, which may vary substantially over the pulsation cycle.

Intriguingly, we found a strong degree of spatial alignment between the high-probability extratidal stars and the leading arm of the simulated tidal stream both in the near- and far-field around the cluster. Unfortunately, there are no SDSS plates covering the in-between regions and so the existence of a stream remains uncertain.

Navin et al. (2016) also found potential extratidal associates for M 3. Following the same approach and reasoning as in Sect. 4.1.1, we exclude previous membership for all eight with high confidence. Two stars can even be discarded based on their highly significant Gaia parallaxes which render them much closer than the cluster (7 kpc and 5.7 kpc compared to 10.2 kpc distance to M 3). The same graphical representation shown for M 13 in Fig. 8 is presented in Fig. 11 for M 3.

Fig. 11. Same as Fig. 8, but for M 3.

4.1.4. M 2 (NGC 7089)

M 2 is a halo (R_GC ∼ 10.1 kpc) GC with a complex chemical-enrichment history. Yong et al. (2014) and Lardo et al. (2016a), for example, suggested an accretion origin of this cluster, which may be the stripped core of a dwarf spheroidal galaxy based on observed anomalous spreads in iron and neutron-capture elements. Kuzma et al. (2016) detected a photometric stellar envelope surrounding M 2 beyond its tidal radius.

Our formalism yielded 77 associates including 11 extratidal candidates for M 2; three of which are uncertain to be truly unbound because of their small projected separations from the GC (< 1.2 r_t). While ten of the stars might be associated from a photometric point of view (see Fig. 12), one candidate is consistently reported by both Gaia and SDSS to be redder than the cluster population at comparable brightness on the MSTO. A further two to three stars reside in the BSS region; at this point we can neither confirm nor reject their BSS nature. Owing to the sparse coverage of SDSS plates around M 2, we cannot draw firm conclusions with respect to the spatial distribution of extratidal halo stars for this GC. However, in the immediate surroundings, we can chemodynamically confirm the finding of Kuzma et al. (2016), who reported on a diffuse stellar envelope that extents out to at least 5 r_t.

Fig. 12. Same as Fig. 7, but for M 2. The CN-strong star from Paper I is labeled “CN”. The red outlier in the upper left panel is outside of the plotting range in the upper right panel as showing it (g0 = 20.71 mag, (g − r)0 = 1.05 mag) would strongly distort the diagram.

The only identified extratidal giant coincides with the parameter range of – and thus was analyzed in – Paper I; it was found to be CN-enhanced, again suggesting a large fraction of 2P stars in the group of escapees (see Sect. 4.5 for further discussion).

4.1.5. M 15 (NGC 7078)

At a slightly lower [Fe/H] than M 92, the cluster M 15 is the most metal-poor cluster for which we could identify associations. The region around M 15 is barely covered by SDSS plates (see Fig. 13). Nevertheless, we find five halo field stars that are a chemodynamical match to this GC. All of them inherit an excellent photometric consistency with the cluster CMD on the MSTO and the blue and red HB. Thus, a cluster origin appears highly feasible.

Fig. 13. Same as Fig. 7, but for M 15.

The star falling in the instability strip is listed as variable in Gaia DR2. Nonetheless, it is not tabulated in either of the works by Drake et al. (2013a,b, 2014), nor in Abbas et al. (2014). We therefore performed our own light curve analysis (Appendix A) and found this star to be a double-mode RR Lyrae of type d. A distance of 9.5 ± 2.3 kpc was deduced, which places the star at a distance in good agreement with that of M 15.

4.1.6. M 53 (NGC 5024) and NGC 5053

M 53 and NGC 5053 are two GCs at a low angular separation of 0.96 deg, even though their real spatial separation is about 1 kpc (Jordi & Grebel 2010). Law & Majewski (2010a) suggested that both clusters could be possible associates of the Sagittarius stream, while Forbes & Bridges (2010) speculate over the possibility of one or both of them being the nucleus of a disrupted dwarf galaxy. The former hypothesis has been refuted with high confidence by Sohn et al. (2018) using Hubble Space Telescope proper motions. Chun et al. (2010) report on extended overdensities for both clusters and claim the detection of a tidal bridge between them. However, Jordi & Grebel (2010) could not reproduce that finding.

Among the extratidal field star population around M 53, we found seven stars that – based on their kinematics and metallicity – may have originated in M 53. As can be seen in Fig. 14, photometrically, two of the stars can be readily associated with the RGB. One star is fainter than the HB in G while being a perfect match in g₀, thus indicating photometric variability. Indeed, Gaia DR2 lists this star as variable and Abbas et al. (2014) classify it as an RR Lyrae star of type c and provide a distance of 24.07 kpc. This appears marginally consistent with the cluster’s distance of 17.9 kpc. Nevertheless, as for the RR Lyrae in Sect. 4.1.3, a comparison is prohibited in light of missing error margins. Our own analysis showed a distance of 20.1 ± 4.8 kpc, which is in agreement with that of M 53 within one error margin. Kundu et al. (2019) also reported on five extratidal RR Lyrae for M 53. Again, due to their restricting themselves to the immediate cluster vicinity, our target is not part of their list.

Fig. 14. Same as Fig. 10, but for M 53.

The remaining four extratidal associations with M 53 cannot be photometrically matched to either the RGB or the HB. Nevertheless, they cover the same colors and magnitudes as the comparably large number of intratidal stars of the BSS population that we attributed from Gaia kinematic properties alone (cyan dots in Fig. 14).

For NGC 5053, we report on 24 plausible extratidal detections, 12 of which have their chemodynamical membership bolstered by a reasonable photometric match to the populations of intratidal members on the blue HB (five stars) and RGB (seven stars, see Fig. 15). Two stars are clearly too red to be attributed to the cluster, while ten others lie in the BSS region of the CMDs with their affiliation nature remaining spurious.

Fig. 15. Same as Fig. 7, but for NGC 5053.

We emphasize that none of the associations are shared between the two clusters considered in this section. This does not necessarily exclude the possibility of a tidal bridge for two reasons: Firstly, we do not model such tidal interactions when simulating the streams. Hence, in the near-field of the clusters, the streams are spatially almost parallel and do not cross. The differences in v_r between the clusters exclude mutual associations. The second, more important reason is that the intercluster field was only covered by a few (∼10) SDSS fibers, thereby tremendously reducing the chance of finding a potential bridge associate.

4.1.7. NGC 4147

NGC 4147 is an outer-halo (R_GC = 21.4 kpc) GC that has been suggested to be associated with the Sagittarius stream by Mackey & van den Bergh (2005). However, Villanova et al. (2016) cast some doubt on this by showing that the cluster is a closer match to the Galactic halo in terms of their respective chemical abundances. More recently, proper motion studies firmly rejected an association since the cluster is on a counter-rotating orbit with respect to Sagittarius (Sohn et al. 2018; Riley & Strigari 2020).

Owing to the faint nature of NGC 4147 (horizontal branch magnitude, V_HB = 17.02 mag, Harris 1996), our approach can merely test targets that are as bright as the RGB of this GC. Of 19 extratidal candidates, we identified three associations that fall on top of the HB of the cluster. The reddest of these three is an RRab star at a distance of 22.0 ± 5.2 kpc⁷. The latter distance renders it consistent with NGC 4147 at 19.3 kpc. Most of the remaining associations (8 stars) were found to be on the upper RGB, whereas two stars below the HB, two stars redward of the RGB, and the five faintest candidates are likely to be false-positive detections judging from their photometric discrepancy from the intratidal Gaia sources (see Fig. 16).

Fig. 16. Same as Fig. 7, but for NGC 4147.

4.2. Associating CN-strong stars with clusters and major merger events

We explored chemodynamical links of the sample of 112 CN-strong field stars from Paper I with the Galactic GC population through method II. Similar efforts have recently been made by Savino & Posti (2019) who used the similar but smaller sample of Martell et al. (2011). We emphasize that only 27 of the 63 cluster-star pairs reported by Savino & Posti (2019) can be directly compared to our study because the remaining ones lie outside of our stellar metallicity restriction of −1.8 <  [Fe/H] <  −1.3 dex⁸ (Paper I). Of these 27, only one is marginally attributable to their sample of best association candidates (< 65% rejection confidence). The reason for the low overall overlap is that we find generally lower probabilities because we use the full distribution of actions instead of marginalized distributions for each of the three components.

Comparing each CN-strong halo star i with each cluster j yields a matrix p_ij with 16 800 entries. Figure 17 depicts this matrix in a representation where both rows and columns are sorted by increasing metallicity. Once again, we highlight the fact that a value for p close to unity does not necessarily imply high confidence in association, but rather that association cannot be rigidly excluded. In case of considerable ambiguities for the orbital parameters of a particular star – for example due to a highly uncertain distance – almost no cluster can be kinematically excluded from being a former host. In such cases, p is entirely driven by [Fe/H] alone. In Fig. 17 those manifest in rows essentially showing broad normal distributions with peak positions at the intersect of the stellar and cluster [Fe/H]. The latter are only occasionally interrupted by low- and high-p columns representing GCs in the outermost regions of the Galaxy and the GCs E 1 and Pal 4, respectively. The latter have loosely constrained orbits themselves, and therefore we do not consider these two clusters in the following.

Fig. 17. Graphical representation of the matrix pij with clusters on the abscissa and CN-strong stars on the ordinate. Both coordinates are sorted by increasing [Fe/H] as indicated by the top and right panels. The five pairs with highest modified confidence, $p_{ij}^{\prime }$ (Eq. (9)), are highlighted by colored circles and presented in the legend.

Inverting the above reasoning implies that pairs (i, j) are good association candidates if they have a high attributed p, while at the same time p_ij/∑_jp_ij is strongly peaked. Hence, each cluster can in principle have an arbitrary number of stellar associations whilst to be considered as part of a good pair, each star should have a strongly limited number of attributable clusters. This leads us to the modified quantity:

$$\begin{array}{c}\hfill p_{\mathit{ij}}^{\prime }=\frac{p_{\mathit{ij}}^2}{{\sum }_j{p}_{\mathit{ij}}}·\end{array}$$(9)

In total, 145 pairs were found to satisfy $p_{ij}^{\prime }>0.05$, whereas only 15 posses values above 0.32. The five pairs with highest p′ are highlighted in Fig. 17 and the whole list can be found in Table C.2. We recover all four associations of CN-strong stars to the GCs M 13 and M 2 that were presented in Sect. 4.1. In case of the association with M 2, there are two additional clusters that come into question given their action integrals (M 22 and ESO 280−06). Given the spatial coincidence with M 2 – which did not enter the present analysis – we favor association with this cluster.

Of the 145 reported star-cluster association pairs, 26 involve GCs that were proposed to have possibly been accreted as part of either Gaia-Enceladus (Myeong et al. 2018) or the Sequoia merger event (Myeong et al. 2019); an observation that represents tentative evidence that the involved CN-strong stars were donated by those galaxies. The pairs are indicated in Table C.2. Moreover, among the 15 strongest ($p{\prime }_{\mathit{ij}}\ge 0.32$), three pairs involve the bona-fide Enceladus clusters M 75 and NGC 1261. Overall, we found that most, namely 23, pairs are accounted for by Enceladus clusters such as M 2, while only three pairs with low association probability (≤0.08) involve Sequoia GCs. We caution that several stars have associations not only with clusters from the two merger events but with other clusters, too. One star, Gaia DR2 603202356856230272, at the same time shows associations with two Enceladus GCs (M 79 and NGC 1851) and one Sequoia cluster (NGC 3201), although with generally low probabilities.

We confirm the finding by Savino & Posti (2019) that a substantial fraction (here 38%) of the CN-strong stars from the sample seemingly cannot be chemodynamically associated with any of the known GCs. We recall here the three reasons these latter authors proposed for this observation under the restrictive assumption that GCs are the only birth site of these chemically peculiar stars: The first, rather unlikely option is that the hitherto unassociated stars could originate from yet-to-be-discovered clusters that are concealed by the high-extinction regions of the Galactic disk or bulge. Secondly, the stars in question could stem from already entirely dissolved clusters. Finally, the stars could have originated from one of the known GCs without retaining the orbital characteristics of the cluster. For this latter case, Savino & Posti (2019) offer two highly plausible explanations involving high ejection velocities due to three-body interactions or an early escape from the cluster followed by a drastic change of the Galactic potential due to, for example, a major merger. In order to address the former scenarios, it is of utmost importance to establish a ground truth by rigidly confirming the CN-strong nature of the targets using spectroscopic follow-up analyses that could reveal other light-element anomalies.

4.3. Associations around CN-strong field stars

Of the 112 investigated stars, method III revealed 81 groups with at least two – possibly, but not necessarily CN-normal – additionally associated stars. Of these groups, 69 comprise five or more stars, whereas 51 still consist of at least ten stars in excess of their respective CN-strong candidate. Naturally, we recover a large number of intratidal targets for the four extratidal CN-strong stars in the vicinity of M 13 and M 2 that were reported in Sect. 4.1.

All remaining CN-strong stars have no obvious direct connection to any cluster. Nevertheless, many of the chemodynamically linked groups of stars can be photometrically attributed to the same stellar population even in the absence of a cluster (see Fig. 18 for an example). Since our selection function is not only spatially inhomogeneous but also strongly biased towards the brighter evolutionary stages, many associated stars coincide with the HB in the CMDs of the latter populations. For candidate variable stars that fall in the instability strip, we performed the analysis outlined in Appendix A and recomputed the Mahalanobis distance using updated values for v_r. This way, we were able to associate 31 already known RR Lyrae stars (Drake et al. 2013a,b, 2014; Abbas et al. 2014) and add one more RRc-type pulsator.

Fig. 18. Graphical representation of the population that was chemodynamically associated with the CN-strong star Gaia DR2 2507434583516170240 (Sgr CNs 3). The color coding and circle sizes have the same meaning as in Fig. 7. Additional red stars indicate the identified RR Lyrae variables (see main text). Left panel: position on the sky on top of rejected stars from the SDSS footprint (gray dots). Middle panel: photometric associations in the Gaia CMD. To guide the eye, isochrones that represent ages from 4 to 13 Gyr (white to dark blue) at the mean [Fe/H] of the population are shown for reference. The CN-strong star is labeled “CN“. Right panel: RR Lyrae G magnitudes and their inferred distances (red stars) in comparison to the spectrophotometric distance of the CN-strong star (black circle, see Sect. 3.2).

Unfortunately, a number of factors prevented a straightforward automated classification for the degree of belief in true photometric association. Among these are the occasionally high rate of obviously spurious associations and the ambiguity of the distance modulus in the absence of significant parallaxes and/or attributed RR Lyrae stars. We present those 17 groups that passed a visual inspection. Their CMDs and spatial distributions can be found in Figs. 18 and B.1 through B.16, while the relevant information about individual stars is presented in Table C.3. Among all CN-strong stars, only Gaia DR2 615481011223972736 and Gaia DR2 634777096694507776 (cf., Figs. B.13 and B.15) were pairwise attributed to each other using method III, which provides evidence that they may share the same birth place.

For five of the groups of associated stars, the CN-strong star was attributed to at least one of the Enceladus clusters through method II. However, four of those only have very low association probabilities, which is why we do not propose an Enceladus membership. The remaining star, Gaia DR2 3696548235634359936, generates the sixth strongest among all star-cluster associations ($p{\prime }_{\mathit{ij}}=0.54$), the involved cluster being NGC 1261. Nonetheless, it is noteworthy that Myeong et al. (2018) classified this cluster as merely a possible member of Enceladus whilst eight additional GCs were labeled probable members.

4.4. Associations to the Sagittarius stream and M 54

We found four of our high-confidence groups of associations from the previous section to coincide with the Sagittarius (Sgr) stream. For further reference, we label the groups Sgr CNs 1 to 4. Fortunately, all four have at least one attributed RR Lyrae star, thus enabling a meaningful 6D phase-space characterization. Figure 18 exemplarily illustrates the tight photometric consistency for Sgr CNs 3, while in Fig. 19 we present various population representations on top of the N-body simulation of the tidal disruption of the Sgr dwarf spheroidal galaxy in a triaxial MW halo by Law & Majewski (2010b). The phase-space comparison places our agglomerations right on top of the dynamically young portion of the leading and trailing arms of the stream. On the basis of the simulated model parameters – under the assumption of true association – the stars became unbound from the Sgr main body no longer than 3 Gyr ago.

Fig. 19. Multi-parametric representation of four populations of chemodynamical associations with CN-strong giants (salmon, yellow, dark red, and light blue circles) on top of the debris model of the Sgr stream in a triaxial halo by Law & Majewski (2010b). The latter is depicted by dots where increasing color lightness denotes a longer elapsed time since becoming unbound from the Sgr main body. Top: 3D position in Cartesian, Galactocentric coordinates. Red-, black-, and violet-dashed lines indicate the positions of the Sun, the Galactic center, and the Sgr main body, respectively. Bottom panels: projections in all six dimensions of positional and kinematic observables. For the bottom panel, we assumed all stars to reside at exactly the same D as inferred from the mean RR Lyrae distance for each population.

From a chemical point of view, the metallicities of the four involved CN-strong stars cover a narrow range from −1.51 to −1.40 dex and thus – within their errors – are fully consistent with the mean [Fe/H] of M 54 (−1.45 dex, e.g., Bellazzini et al. 2008). This observation suggests that four chemically altered giants and potentially a considerable fraction of their associates identified here originated from the massive GC M 54. Indeed, one of the four involved CN-strong stars – namely Gaia DR2 2507434583516170240, giving rise to the association Sgr CNs 3 – has already been linked to this cluster based on its position in action space (method II), though with very low confidence (0.08). We suspect that the reason for the latter is to be found in the main flaw of the assumptions in Sect. 4.2, that is, the fact that actions of potential cluster escapees do not necessarily have to be exactly identical to the cluster itself; otherwise, they would never have been able to become unbound in the first place.

4.5. The fraction of chemically altered stars amongst bona-fide escapees

In Paper I we investigated the fraction of CN-strong stars in our sample of analyzed halo stars to estimate the fractional GC contribution to the stellar inventory of the Galactic halo, $f_{\mathrm{h}}^{\mathrm{GC}}$. The outcome of this assessment, among other factors, depends on the assumed fractions of chemically normal 1P stars with respect to the total number of escapees (N_1P/N_tot), resulting in a corridor for the overall GC contribution to the halo between 8.5% and 12.5%. As discussed in detail in Paper I, several values for N_1P/N_tot have been proposed for the still bound stellar populations of GCs. For instance, Carretta et al. (2009) and Bastian & Lardo (2015) reported 50% and 32%, respectively, irrespective of GC properties. Milone et al. (2017) found a dependency on the present-day mass of the clusters. Here, for the first time, we elaborate on the average N_1P/N_tot in the collection of escaped stars in the Galactic halo by employing our established chemodynamical links.

Methods I and III revealed new potential GC escapees in the halo field star population that are either directly chemodynamically linked to a cluster or – adopting the assumption that GCs are the only production sites of CN-strong stars – can be associated with a CN-strong star. Here, we account for the fact that the candidates cover a broad range of evolutionary stages, while Paper I focused exclusively on the RGB. Hence, the newly discovered sample was cleaned to match the selection criteria outlined in Sect. 2.1 of Paper I. Briefly, the targets were split into metallicity bins and the exact same fiducial regions for the RGB in a log g–(g−r)₀ diagram were employed. Most RGB stars with an established chemodynamic link in this study are lost for the present analysis as they fall outside of the stringent metallicity corridor applied in Paper I (−1.8 <  [Fe/H] <  −1.3 dex⁹). As a result, our calculations exclude, for example, the metal-poor cluster M 92 with its many associations on the RGB. Despite reducing the statistical significance, restricting the sample to the RGB confers the coincidental advantage that spurious associations can be considered less likely as the involved stars are brighter than, for example, the upper MS. This results in higher precision in the astrometric quantities [Fe/H] and v_r and therefore stronger discriminatory power.

Figures 20 and 21 illustrate the subsample obtained in this way for methods I and III. We note that we see no reason to suspect that 1P and 2P stars would have been subject to different selection functions in the SDSS fiber placement, which was solely based on pre-selection from SDSS photometric broad-band filters which in turn are not sensitive to light-element anomalies. Given that method III cannot firmly distinguish between former 1P GC stars and stars stripped from the Sgr dwarf spheroidal galaxy, we excluded the four involved CN-strong stars and their associations. The same argument holds for the population that was tentatively associated with Gaia-Enceladus (Sect. 4.3), which is why we exclude the CN-strong star and its one associated giant obeying the selection criteria, too. Having identified four bona-fide 1P and four 2P stars through method I and 48 1P and 12 2P stars from method III, we calculated N_1P/N_tot = 50.0 ± 16.7% and $80.2_{-5.2}^{+4.9}\%$, respectively. The values and error margins are taken from the median, 15.9th, and 84.1th percentiles of distributions that were generated in a Monte Carlo (MC) analysis using 10⁷ random draws from Poisson distributions for the involved counting statistics. The two resulting distributions can be found in Fig. 22 where we also show the findings reported for bound cluster stars that were mentioned earlier in this section.

Fig. 20. Gaia CMD (left) and Kiel diagram (right) of all association candidates from method I. For the absolute magnitude, MG, we computed a distance modulus under the assumption that all stars reside at exactly the same distance as their potential former hosts. Dark gray to light gray colors indicate the stars’ metallicity according to the color bar, while the point size mirrors the calculated association probability. The metallicity-dependent polygons used for pre-selection in Paper I are depicted in the right panel. Associated stars which at the same time obey these selection cuts are shown in blue. Red circles depict the four CN-strong stars that are likely associated with M 13 and M 2, respectively.

Fig. 21. Same as Fig. 20 but for associations around CN-strong stars (method III). Here, MG was computed using the spectrophotometric distances of the CN-strong stars as outlined in Sect. 3.2. Orange and green circles show the associations with the Sgr stream and Gaia-Enceladus that obey the selection criteria of Paper I but were not considered when estimating the fraction N1P/Ntot.

Fig. 22. Result of an MC simulation with 107 realizations for the deduced fraction of first-generation stars among all potential GC escapees obeying the selection function in Paper I (Figs. 20 and 21). Here, the blue and red histograms represent the findings from methods I and III, respectively. For comparison, we indicate the fractions for bound cluster stars as reported by Carretta et al. (2009, solid vertical line), Bastian & Lardo (2015, dashed vertical line), and the lowest value from Milone et al. (2017, dash-dotted vertical line).

Interestingly, in the near-field around clusters (i.e., using method I), the first-generation fraction appears to balance the second-generation fraction, and is therefore broadly consistent with the present-day cluster findings by Carretta et al. (2009) and Bastian & Lardo (2015), and with the majority of the analyzed clusters of Milone et al. (2017). Moving away from the GCs, that is, to regions in the halo without a direct spatial link to any GC, N_1P/N_tot seems to increase considerably. While it is tempting to claim a solid finding, in particular in light of the low attributed statistical errors from method III, we discuss here two effects that may lead to biases in the estimates.

The first possibility is that some of the CN-strong stars could have originated in a cluster that was part of an accreted, now entirely disrupted dwarf galaxy (other than Sgr, Gaia-Enceladus, and Sequoia) with its own system of GCs. Therefore, much like already seen in the five CN-strong stars of the Sgr stream and (possibly) Gaia-Enceladus, the stars associated with those CN-strong stars would not only share the chemodynamical information of the cluster escapees but that of the stars of the galaxy as well¹⁰. Vice versa, while being characterizable as ex situ, associates of the CN-strong stars would not necessarily be of GC origin. As a consequence, the number of chemically normal stars that are erroneously classified as 1P stars would be artificially increased, which in turn would lead to a higher N_1P/N_tot. Evidence for CN-bimodalities amongst the dSph field star populations is sparse (e.g., Smith & Dopita 1983; Lardo et al. 2016b; Norris et al. 2017; Salgado et al. 2019). However, several GCs in the MW have been associated with dwarf-accretion events. It is equally likely that GCs that were part of dSphs have dissolved into their field star populations (prior to their merger with the MW; see Malhan et al. 2019) meaning that the presence of CN-strong stars within dSphs could be expected, provided that the GCs of dSphs follow the same trends as the in situ population of the MW (e.g., Hendricks et al. 2016).

Another source of increased false-positive rates for 1P stars are spurious associations from the in situ population of halo stars. However, this latter possibility seems rather unlikely given the independent finding of tight photometric consistency even outside of the selection bounds, whereas halo interlopers could in principle occupy any position in the CMD.

Bearing in mind all the caveats mentioned above, taken at face value, our high N_1P/N_tot for the halo far away from GCs could imply that 1P stars were preferentially lost during early cluster dissolution while 2P stars were more easily retained. This would provide support for early mass-loss scenarios that commonly assume that only 1P stars are lost (e.g., D’Ercole et al. 2008, see also the review by Bastian & Lardo 2018 and the discussions in Paper I). Alternatively, assuming that GC mass loss affects both 1P and 2P stars to the same extent (Kruijssen 2015) and considering the inverse mass dependency of the now-observed occurrence rate of 1P stars following Milone et al. (2017), low-mass clusters could have contributed the majority of the former GC stars now found in the halo. For the fraction of halo stars donated by GCs, the findings presented here, in concert with the formalism outlined in Sect. 5.1. of Paper I, yield $f_{\mathrm{h}}^{\mathrm{GC}}=10.5\pm 0.7\%$ and 11.8 ± 0.2%, assuming either the 1P fraction for the cluster vicinities or for the general halo, where we adopted an early mass-loss rate of 56% and an average present-day cluster mass of 2.2 × 10⁵ M_⊙ (see Paper I, for detailed discussions).

5. Summary and conclusions

Here, we investigate the connection between GCs and the stellar component of the Galactic halo in order to contribute to our understanding of how this old component of the Galaxy was formed. To this end, we explore chemodynamical links between halo field stars and the present-day cluster population using combined data from the SDSS/SEGUE and BOSS surveys as well as from the Gaia mission. We created a dataset that allows for the characterization of about 3 × 10⁵ stars through their full phase-space vectors and chemistry (by means of metallicity). As realistic error budgets are key to our probabilistic approaches, uncertainties on SDSS radial velocities and metallicities were reassessed using well-constrained cluster populations and as a consequence considerably increased. The most important ingredient of our study is the characterization of 112 giant stars in the halo that in our previous study (Koch et al. 2019) were found to be chemically peculiar in the sense that they show atypically high nitrogen abundances. In the present investigation, we build upon the assumption that the environmental conditions required for the production of such peculiar chemical characteristics are unique to GCs, the only sites where an abundant occurrence of these 2P stars has been found so far.

Our analysis focuses on three techniques that can be distinguished by the stringency they impose on the precision of the input data. The first approach is restricted to the immediate volumes around current cluster positions, where a solid distance estimate – made using significant parallaxes – is not strictly required. In this case, field stars were associated using only six parameters (two spatial and three velocity coordinates as well as [Fe/H]). The confidence in true association was independently solidified by the fact that the associated stars show a remarkably close match to the still bound populations of stars. The second method uses action integrals to tie CN-strong stars to clusters under the presumption that cluster escapees retain their kinematic memory. A spatial coincidence with the clusters is not required at the expense of a much stronger demand in terms of well-constrained distance estimates. This prohibits a wider applicability for distant halo stars in excess of the 112 CN-strong stars, for which we inferred spectrophotometric distances. Nevertheless, we discuss several persistent systematic uncertainties that may have considerable impact on such measurements. The third and final tagging approach adopted throughout this work was used to single out stellar populations in the near-field of CN-strong stars that share similar chemodynamical properties. It may therefore be feasible that those populations stem from the same birth place, that is, GCs.

In total, we established direct chemodynamical links in the near field of eight clusters (NGC 4147, M 53, M 3, M 13, M 92, M 15, M 2) for 789 stars. While the projected distances of 638 stars are found to reside within the clusters’ tidal radii – thereby rendering them likely bound cluster members – 151 can be denoted extratidal. Among the escapees, we qualitatively recover structures that resemble the results of previous photometric studies identifying extratidal envelopes (e.g., Jordi & Grebel 2010). The associations are spread over a wealth of evolutionary stages and metallicities (e.g., four of the eight GCs are more metal-poor than −2 dex). Nonetheless, eight giant stars share the same selection function as required in Koch et al. (2019), and four of these are recognized as being CN-strong. From this, based on the low number statistics at hand, we estimate the fraction of 1P stars among the close-by escapees to be N_1P/N_tot = 50.0 ± 16.7%, which is broadly in line with the results of studies that were dedicated to estimating this fraction in nowadays observable cluster populations (Carretta et al. 2009; Bastian & Lardo 2015; Milone et al. 2017).

We further report on 145 possible (p >  0.05) pairs of CN-strong stars in the halo field star population and individual GCs. Of these, 15 are presented as probable star-to-cluster connections (p ≥ 0.32). About 18% of the association pairs involve clusters that were proposed to be attributable to the major merger events Gaia-Enceladus and Sequoia. Three of the high-probability pairs involve Enceladus GCs (M 75 and NGC 1261). We show that 38% of the CN-strong stars cannot be linked to any surviving cluster, which is in marginal agreement with the finding of ∼50% by Savino & Posti (2019). Among the possible explanations is a scenario whereby these unassociated stars originate from already destroyed clusters.

The third tagging approach revealed 17 populations of stars that share the same portion of the chemodynamical information space as CN-strong stars. For nine such agglomerations of stars, we are able to attribute already known and newly discovered RR Lyrae stars, for which our refined approach for the distance inference provides more reliable distance estimates than the spectrophotometric distances with their attached caveats. Of particular interest, using this method, four CN-strong targets together with their associated stars are found to be probable members of the Sgr stream. Based on the metallicity overlap with M 54, a GC in the central part of the Sgr dwarf galaxy, we suggest that the four stars may have been stripped from this cluster.

From the identified RGB stars among the stellar populations around CN-strong stars without spatial coincidence with GCs, we again determined N_1P/N_tot and find a value of $80.2_{-5.2}^{+4.9}\%$, which should be regarded as tentative because the statistical uncertainties do not account for all the potential caveats discussed. Nevertheless, it is tempting to argue that a substantial increase of N_1P/N_tot from the vicinity of clusters to the overall halo field points towards either a preferential loss of 1P stars during early cluster dissolution or the birth of a considerable fraction of former GC stars of the halo in since-dissolved clusters. In any case, our findings provide a powerful observational benchmark for theoretical studies of cluster disruption.

As already emphasized in our earlier work, unambiguous confirmation that the identified CN-strong stars are indeed agents of the GC populations showing light-element anomalies remains to be attained. Evidence for imprints from high-temperature proton burning could be gathered through spectroscopic analyses of other key tracer elements such as O, Na, Mg, and Al. Another layer of complexity is added by the open question of whether GCs are genuinely the only sites providing the necessary conditions for the development of abundance anomalies. Bekki (2019), for example, presented a scenario that could explain the substantial fraction of N-rich stars found amongst the Galactic bulge field stars; namely, that these N-rich stars may have formed in the field from an ISM that has been pre-enriched by AGB ejecta. This latter author further envisions a similar mechanism for the halo. However, our – admittedly tentative – observation of an overabundance of CN-strong stars in the cluster surroundings as compared to the halo field challenges this explanation and suggests a mixture of both GC contribution and in situ formation in the case of the halo.

Our method that is capable of identifying peculiar stars (Paper I) was restricted to the metallicity range −1.8 <  [Fe/H] <  −1.3 dex. Nevertheless, most of the associates reported in the present study, in particular the ones directly surrounding clusters, are more metal-poor than the lower [Fe/H] limit used in Paper I. For instance, the very metal-poor GC M 92 with its many RGB associates both in the neighboring extratidal envelope and at larger separations remained untouched by our previous approach. Hence, a high-resolution follow up of these moderately faint (14.5 <  G <  17.5 mag) targets is strongly called for in order to determine the extratidal N_1P/N_tot in the metal-poor regime. Such investigations are paramount to observationally test simulations predicting that the fraction of 1P stars in the halo should increase with decreasing metallicity (Reina-Campos et al. 2020). In addition, not being restricted by the accessibility of CN band indices, high-resolution follow-up studies could target our associated MS, MSTO, SGB, and HB targets as well. This would considerably increase the statistical basis of the fractional estimates. The availability of improved number statistics would allow trends with other parameters such as the cluster age to be investigated (e.g., Martocchia et al. 2019). Finally, in the very near future, ongoing and upcoming large-scale spectroscopic surveys such as APOGEE (Majewski et al. 2017), WEAVE (Dalton et al. 2012), and 4MOST (de Jong et al. 2012; Helmi et al. 2019; Christlieb et al. 2019) will substantially enlarge the sample of chemodynamically characterized stars in the Galactic halo. The exploration of those vast datasets that are not restricted to a handful of chemical elements will ultimately lead to a much refined picture of the connection of GC stellar populations and the Milky Way halo.

Acknowledgments

M.H., A.K., and E.K.G. gratefully acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subprojects A03, A05, A08, and A11). Z.P. acknowledges the support of the Hector Fellow Academy. The authors are grateful to G. Parmentier for fruitful discussions. We highly appreciate the constructive report by the anonymous referee. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182 and AST-1313422.

References

Appendix A: RR Lyrae analysis

RR Lyrae stars are pulsating variables residing on the horizontal giant branch within the instability strip. Their pulsation periods are tightly connected to their luminosity through period–metallicity–luminosity relations (see, e.g., Catelan et al. 2004; Muraveva et al. 2018), which makes them invaluable distance indicators within the Local Group. Furthermore, they are part of the old stellar population with ages > 10 Gyr due to their exclusive occurrence in old stellar systems such as globular clusters (VandenBerg et al. 2013). RR Lyrae variables are divided into three classes based on the pulsation mode; fundamental mode stars (RRab), first overtone mode pulsators (RRc), and double-mode variables (RRd, pulsating simultaneously in the fundamental and first overtone).

A few hundred stars in our analysis fall inside the fundamental mode RR Lyraes’ instability strip boundaries provided by Clementini et al. (2019), intrinsically making them candidates for possible variability. In order to investigate their periodic alternation, we retrieved their photometry from the time domain Catalina sky survey (CSS, Drake et al. 2009). In addition, we used datasets assembled using the data from the CSS (Drake et al. 2013a,b, 2014; Abbas et al. 2014) as a bona fide catalog of RR Lyrae stars in the Galactic halo. The crossmatch of the RR Lyrae candidates with the halo RR Lyrae sample resulted in matches for 65 cases with determined pulsation periods. To verify the nonvariability of the remaining stars, we elaborated on their possible periodicity using the upsilon package¹¹ (Kim & Bailer-Jones 2016), which searches for periodic behavior among photometric data and provides a classification and light curve parameters of identified variables. To ensure reliable classification, we required at least 50 % class probability for stars identified as RR Lyrae pulsators which were fulfilled by nine stars that we denote as newly discovered RR Lyrae variables (three RRab, five RRc and one RRd).

The one RRd variable in our sample is a newly discovered double-mode RR Lyrae star with a dominant first overtone mode (pulsation period P_1O = 0.4040992 d) and a secondary fundamental mode at P_F = 0.5417994 d. The period ratio between both modes is P_1O/P_F = 0.7458465, which together with the long fundamental mode period suggests a low metallicity between −2.0 and −2.5 dex (see Fig. A.1 for details) in line with M 15’s [Fe/H].

Fig. A.1. Period ratio vs. fundamental pulsation period (errors on depicted quantities are smaller than the point size) for a discovered RR Lyrae star (gray point) compared to the pulsation models (blue and red colors, Smolec & Moskalik 2008) calculated for a fixed mass (0.7 M⊙), metallicity (−2.0 and −2.5 dex), and a grid of stellar parameters (L and Teff) typical for RR Lyrae stars.

The assembled RR Lyrae light curves (both known and newly discovered) were phased using the stars’ ephemerids and decomposed using the Fourier series to determine the time of maximum brightness, the mean magnitude, and the pulsation amplitude:

$$\begin{array}{c}\hfill m(t)=A_0+{\displaystyle \sum _{k=1}^n{A}_k·\cos (2\pi k\vartheta +{\phi }_k),}\end{array}$$(A.1)

where ϑ represents the phase function defined as (MJD − M₀)/P, MJD denotes the observation in the Modified Julian Date, and P and M0 are ephemerids of a given star (pulsation period and time of brightness maximum). A_k and φ_k stand for amplitudes and phases, respectively, with n denoting the degree of the Fourier series, and A₀ representing the mean magnitude. Each phased light curve was visually inspected to verify the ephemerids and pulsation amplitudes derived from the Fourier fits.

The absolute magnitudes of RR Lyrae stars in the visual part of the spectrum strongly correlate with their metallicities (see, e.g., Catelan et al. 2004; Muraveva et al. 2018), thus acquiring homogenous absolute magnitudes for all three RR Lyrae classes in our sample can be troubling due to the lack of reliable metallicities (cf., e.g., Smolinski et al. 2011; Hanke et al. 2018; Fabrizio et al. 2019). Hence, for the distance estimates we decided to assume a single absolute magnitude for the entire RR Lyrae sample and allowed a large dispersion M_V = 0.6 ± 0.5 mag.

The single value for the absolute magnitude of all RR Lyrae variables in our sample was derived using a sample of RR Lyrae stars from Muraveva et al. (2018), Gaia parallaxes, 3D dust maps from Green et al. (2019), and mean V-band magnitudes from the ASAS survey (Pojmanski 2002). For the aforementioned quantities, we ran an MC error simulation of the distance modulus assuming two different offsets in Gaia parallaxes (−0.057 mas, and −0.029 mas, respectively, Muraveva et al. 2018; Lindegren et al. 2018), and estimated absolute magnitudes for individual RR Lyrae variables. The ⟨M_V⟩ magnitudes clustered around 0.68 ± 0.35 mag and 0.56 ± 0.37 mag for Gaia parallax offsets, −0.057 mas and −0.029 mas, respectively. Pursuing a conservative approach, we used M_V = 0.6 mag with a large dispersion of 0.5 mag to account for the offset in both absolute magnitude values derived from parallaxes.

The distances for our sample of RR Lyrae variables were calculated using the MC simulation with 1000 realizations assuming Gaussian error distributions in the reddening (using dust maps from Schlafly & Finkbeiner 2011), absolute magnitudes, and apparent magnitudes¹². The resulting distances come with errors of around 20−25%, larger than generally reported based on optical data for RR Lyrae stars (around 8%, Neeley et al. 2017). This is mainly caused by the conservative assumption for absolute RR Lyrae magnitudes.

Periodic luminosity changes are driven by the variation of the stellar radius and effective temperature. This inevitably affects the radial velocity measurements that comprise a combination of systemic and pulsation velocity. The pulsation effect can be removed from a single radial velocity measurement using radial velocity templates (Sesar 2012)¹³ that scale with their photometric counterparts.

We note that for RRc and RRd type RR Lyrae stars, Sesar (2012) does not provide a scaling relation between the amplitudes of v_r and photometric light curves. Thus, we used a sample collected by Sneden et al. (2017) for radial velocity amplitudes (based on the Hα line) of the first overtone pulsators in the field and compared them with their photometric amplitudes using data from the ASAS survey (Pojmanski 2002). Applying the latter comparison we obtained the following scaling relation between the radial and photometric amplitudes for the first overtone pulsators:

$$\begin{array}{c}\hfill A_{\mathrm{rv}}^{\mathrm{H}\alpha }=95.8A_V-5.2.\end{array}$$(A.2)

The aforementioned relation was also used for the determination of the systemic velocity of the one RRd star in our sample.

In order to determine the systemic velocity of RR Lyrae stars in our sample we extracted single-epoch spectra from the SDSS and measured the radial velocity for the Hα line by cross-correlation with a synthetic template¹⁴ using the iSpec package routines (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019) and MOOG (Sneden 1973, February 2017 version). Individual radial velocity measurements were phased using the stars’ ephemerids with MJDs centered in the middle of the exposure. We scaled the Hα radial velocity template based on the stars amplitude and subtracted the pulsation velocity from the measured radial velocity. The final systemic velocity for a given RR Lyrae star was estimated through the median of individual, pulsation-corrected radial velocity measurements. Errors on the systemic velocities were estimated following the prescriptions by Sesar (2012).

Appendix B: Additional figures

Fig. B.1. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 4018168336083916800 (Sgr CNs 1).

Fig. B.2. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2910342854781696 (Sgr CNs 2).

Fig. B.3. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2464872110448024960 (Sgr CNs 4).

Fig. B.4. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 1020747223263134336.

Fig. B.5. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 1489389762267899008.

Fig. B.6. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2566271722756610304.

Fig. B.7. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3696548235634359936. We note that this group of stars may be associated with Gaia-Enceladus through its association with NGC 1261 (Sects. 4.2 and 4.3).

Fig. B.8. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3707957627277010560.

Fig. B.9. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3735130545329333120.

Fig. B.10. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3896352897383285120.

Fig. B.11. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3906396833023399936.

Fig. B.12. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 603202356856230272.

Fig. B.13. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 615481011223972736.

Fig. B.14. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 615812479620222592.

Fig. B.15. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 634777096694507776.

Fig. B.16. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 777803497675978368.

Appendix C: Additional tables

All Tables

All Figures

Fig. 1. Precision and accuracy validation of [Fe/H]SSPP for M 13 member stars. Lower left panel: proper-motion-cleaned CMD in the GaiaG, GBP, and GRP passbands. For reference to the lower right panel, GBP − GRP is additionally indicated by colored symbols. Right panels: marginalized distribution of deviations (top) and behavior with evolutionary state (bottom) of SSPP metallicities with respect to the literature value: [Fe/H]M 13 = −1.53 dex. Red points and error bars denote the means and standard deviations in bins of 1 mag width in G, whereas underlying gray bars represent the median of the SSPP uncertainties on [Fe/H].

In the text

Fig. 2. Upper left: same CMD as in Fig. 1, but with color coding indicating the brightness (see lower right color bar); the same color code is used in all other panels. Upper right: spatial distribution of the M 13 validation sample. The cluster core, half-light radius (1.69 arcmin), and tidal radius (21 arcmin) are represented by a black star, red circle, and black-dashed line, respectively. We emphasize that the tidal radius together with ϖ and vr were the only selection criteria applied to investigate proper motion systematic errors (main text). Middle panels: absolute (left) and relative (right) deviation in proper motion of the sample from the cluster mean value. Lower panel: squared Mahalanobis distance, M2, from proper motions only versus projected distance from the cluster center. In analogy to the 1σ, 2σ, and 3σ significances of a normal distribution, dash-dotted, dashed, and dotted lines represent the corresponding p-values (see legend). The red line denotes the scaling relation introduced in Eq. (4).

In the text

Fig. 3. Five-parameter chemodynamical association criteria used to constrain the association with the GC M 13. The rejected back- and foreground population of stars with SSPP parameters and Gaia DR2 motions is indicated with gray density contours. Brown dots indicate the simulated stream (see main text). Intratidal and extratidal stars are depicted with blue and colored circles (see also Fig. 7), respectively. Circle sizes indicate the association probability and colors resemble the projected distance to the cluster center (see legend and color bar). The marginalized distributions of each coordinate on the abscissa are shown on top of the two-dimensional representations. The adopted spatial prior as well as the stars that are rejected due to its usage (gray circles with no colored counterpart) are illustrated in the upper right-hand panel.

In the text

Fig. 4. Spatial distribution of GCs (blue squares) and CN-strong halo giants (colored according to their distance) in Galactic coordinates (left) and the Cartesian Galactocentric frame (right). The dashed line in the left panel denotes the celestial equator. The SDSS footprint area is depicted in white while the Gaia all-sky map is illustrated in the background for orientation. Image credit. Gaia Data Processing and Analysis Consortium (DPAC); A. Moitinho/A. F. Silva/M. Barros/C. Barata, University of Lisbon, Portugal; H. Savietto, Fork Research, Portugal.

In the text

Fig. 5. Toomre diagram for our CN-strong giants (colored circles) and the Galactic GC population (blue squares). The latter is only shown in the regime −1.8 <  [Fe/H] <  −1.3 dex, representing the coverage of the analyzed stars. The color coding of the CN-strong stars is the same as in Fig. 4. For better comparability, an additional Gaussian kernel density estimate of the GC sample is provided where increasing color darkness resembles increasing density. To show the clear contrast to the Galactic disk population, we also indicate all targets from our Gaia – SDSS crossmatch that obey −5° < b <  5° and have a spectrophotometric distance as determined by SDSS of less than 1 kpc (gray dots). Moreover, the dashed line encircles the region of absolute velocities of less than 100 km s−1 with respect to the local standard of rest.

In the text

	Fig. 6. Distribution of Monte Carlo samples for the integrals of motion for M 70 (red) and the star Gaia DR2 3833963854548409344 (blue). The marginalized two-dimensional and one-dimensional kernel density estimates are indicated as contours in the frames and histograms at the top, respectively. For illustrative purposes, the histograms were normalized by their maxima.
In the text

Fig. 7. Parameter distributions of intra- (blue) and extratidal (black to light yellow) stars associated with M 13 on top of the field star population covered by the SDSS footprint (gray dots). The sizes of the colored circles are directly proportional to the membership probability P(A|B) as indicated in the legend, while the color of the extratidal candidates denotes the projected distance from the cluster core in units of rt (see color bar). Upper panels: CMDs from Gaia (left) and dereddened SDSS (right) bands. Cyan dots indicate proper motion and ϖ-selected Gaia sources without spectroscopic SDSS counterparts (see main text for details). Typical photometric errors are specified in each case by error bars on the left-hand side. For a subset of bright intratidal stars, reliable color information from SDSS is not available, which is why those targets were omitted in the right CMD. Red lines are 12 Gyr PARSEC isochrones (Marigo et al. 2017) matching the parameters from Table 1 with an applied manual shift in color for the Gaia CMD. The CN-strong giants classified in Paper I are labeled “CN” in both panels. Lower panels: distribution on the sky (left) and a zoomed-in view of the immediate cluster vicinity (right). The black-dashed ellipse indicates the tidal radius, rt, and the direction to the Galactic center is marked by a gray arrow. Properties for the integrated cluster orbit and the simulated stream are indicated by the blue curve and brown dots, respectively.

In the text

	Fig. 8. Comparison of the chemodynamical (left two panels) and spatial (third panel) distributions and the CMD (right) for M 13 between Navin et al. (2016; cyan stars) and the present study (same color and size scheme as in Fig. 7). In the proper motion diagram, the values from the original study are connected by red lines with their Gaia DR2 counterparts (red stars).
In the text

	Fig. 9. Same as Fig. 7, but for M 92.
In the text

	Fig. 10. Same as Fig. 7, but for M 3. The approximate color range of the instability strip from Clementini et al. (2019) is indicated in the upper left panel by vertical dashed lines and the pulsating variable (see main text) is labeled “RRab”. We note that the illustrated color and brightness for the RR Lyrae star are merely mean quantities, which may vary substantially over the pulsation cycle.
In the text

	Fig. 11. Same as Fig. 8, but for M 3.
In the text

	Fig. 12. Same as Fig. 7, but for M 2. The CN-strong star from Paper I is labeled “CN”. The red outlier in the upper left panel is outside of the plotting range in the upper right panel as showing it (g0 = 20.71 mag, (g − r)0 = 1.05 mag) would strongly distort the diagram.
In the text

	Fig. 13. Same as Fig. 7, but for M 15.
In the text

	Fig. 14. Same as Fig. 10, but for M 53.
In the text

	Fig. 15. Same as Fig. 7, but for NGC 5053.
In the text

	Fig. 16. Same as Fig. 7, but for NGC 4147.
In the text

	Fig. 17. Graphical representation of the matrix pij with clusters on the abscissa and CN-strong stars on the ordinate. Both coordinates are sorted by increasing [Fe/H] as indicated by the top and right panels. The five pairs with highest modified confidence, $p_{ij}^{\prime }$ (Eq. (9)), are highlighted by colored circles and presented in the legend.
In the text

Fig. 18. Graphical representation of the population that was chemodynamically associated with the CN-strong star Gaia DR2 2507434583516170240 (Sgr CNs 3). The color coding and circle sizes have the same meaning as in Fig. 7. Additional red stars indicate the identified RR Lyrae variables (see main text). Left panel: position on the sky on top of rejected stars from the SDSS footprint (gray dots). Middle panel: photometric associations in the Gaia CMD. To guide the eye, isochrones that represent ages from 4 to 13 Gyr (white to dark blue) at the mean [Fe/H] of the population are shown for reference. The CN-strong star is labeled “CN“. Right panel: RR Lyrae G magnitudes and their inferred distances (red stars) in comparison to the spectrophotometric distance of the CN-strong star (black circle, see Sect. 3.2).

In the text

Fig. 19. Multi-parametric representation of four populations of chemodynamical associations with CN-strong giants (salmon, yellow, dark red, and light blue circles) on top of the debris model of the Sgr stream in a triaxial halo by Law & Majewski (2010b). The latter is depicted by dots where increasing color lightness denotes a longer elapsed time since becoming unbound from the Sgr main body. Top: 3D position in Cartesian, Galactocentric coordinates. Red-, black-, and violet-dashed lines indicate the positions of the Sun, the Galactic center, and the Sgr main body, respectively. Bottom panels: projections in all six dimensions of positional and kinematic observables. For the bottom panel, we assumed all stars to reside at exactly the same D as inferred from the mean RR Lyrae distance for each population.

In the text

Fig. 20. Gaia CMD (left) and Kiel diagram (right) of all association candidates from method I. For the absolute magnitude, MG, we computed a distance modulus under the assumption that all stars reside at exactly the same distance as their potential former hosts. Dark gray to light gray colors indicate the stars’ metallicity according to the color bar, while the point size mirrors the calculated association probability. The metallicity-dependent polygons used for pre-selection in Paper I are depicted in the right panel. Associated stars which at the same time obey these selection cuts are shown in blue. Red circles depict the four CN-strong stars that are likely associated with M 13 and M 2, respectively.

In the text

	Fig. 21. Same as Fig. 20 but for associations around CN-strong stars (method III). Here, MG was computed using the spectrophotometric distances of the CN-strong stars as outlined in Sect. 3.2. Orange and green circles show the associations with the Sgr stream and Gaia-Enceladus that obey the selection criteria of Paper I but were not considered when estimating the fraction N1P/Ntot.
In the text

Fig. 22. Result of an MC simulation with 107 realizations for the deduced fraction of first-generation stars among all potential GC escapees obeying the selection function in Paper I (Figs. 20 and 21). Here, the blue and red histograms represent the findings from methods I and III, respectively. For comparison, we indicate the fractions for bound cluster stars as reported by Carretta et al. (2009, solid vertical line), Bastian & Lardo (2015, dashed vertical line), and the lowest value from Milone et al. (2017, dash-dotted vertical line).

In the text

	Fig. A.1. Period ratio vs. fundamental pulsation period (errors on depicted quantities are smaller than the point size) for a discovered RR Lyrae star (gray point) compared to the pulsation models (blue and red colors, Smolec & Moskalik 2008) calculated for a fixed mass (0.7 M⊙), metallicity (−2.0 and −2.5 dex), and a grid of stellar parameters (L and Teff) typical for RR Lyrae stars.
In the text

	Fig. B.1. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 4018168336083916800 (Sgr CNs 1).
In the text

	Fig. B.2. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2910342854781696 (Sgr CNs 2).
In the text

	Fig. B.3. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2464872110448024960 (Sgr CNs 4).
In the text

	Fig. B.4. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 1020747223263134336.
In the text

	Fig. B.5. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 1489389762267899008.
In the text

	Fig. B.6. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 2566271722756610304.
In the text

	Fig. B.7. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3696548235634359936. We note that this group of stars may be associated with Gaia-Enceladus through its association with NGC 1261 (Sects. 4.2 and 4.3).
In the text

	Fig. B.8. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3707957627277010560.
In the text

	Fig. B.9. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3735130545329333120.
In the text

	Fig. B.10. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3896352897383285120.
In the text

	Fig. B.11. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 3906396833023399936.
In the text

	Fig. B.12. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 603202356856230272.
In the text

	Fig. B.13. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 615481011223972736.
In the text

	Fig. B.14. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 615812479620222592.
In the text

	Fig. B.15. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 634777096694507776.
In the text

	Fig. B.16. Same as Fig. 18 but for associations to the CN-strong star Gaia DR2 777803497675978368.
In the text