© ESO 2019

1. Introduction

The dynamical evolution of a stellar cluster is strongly affected by internal collisional dynamics, that is, two-body gravitational encounters between stars (see Spitzer 1987, for an in-depth discussion in the case of globular star clusters), stellar-wind mass loss, and external interactions such as encounters with molecular clouds (Spitzer 1958; Gieles et al. 2009) or Galactic disc and spiral arms (Binney & Tremaine 1987). These phenomena modify the cluster structure and internal distribution of the stars along with the cluster lifetime and eventually dissolve the cluster. They are so effective that dynamical models predict that a cluster with the typical mass of a bound open cluster (OC), ≲10⁴ M_⊙, would be dissolved on a timescale of ≲1 Gyr (see e.g. Baumgardt & Makino 2003; Lamers et al. 2005; Bastian et al. 2008).

However, about 20% of the known OCs have ages older than 1 Gyr (Dias et al. 2002; Kharchenko et al. 2013). Although they are a minority, they are key to investigating the dynamical evolution of clusters and the balance of self-gravity, internal dynamics, and external tidal forces. Several features distinguish these old OCs from the bulk of the OCs population. They are systematically more massive than the younger bound clusters. In particular, although these clusters are unambiguously related to the Galactic disc because of their age and chemical composition, they are also often found at high altitude above the Galactic plane, away from the disruptive influence of the disc. Even though, strictly speaking, it would be hard to exclude that some clusters have been acreeted from external galaxies in merging processes.

NGC 2682 (M 67) is one of the nearest old OCs (3.6 Gyr, Bossini et al. 2019) and is located at about 860 pc from the Sun (Cantat-Gaudin et al. 2018a), to which it is quite similar in terms of age and initial chemical composition (e.g. Liu et al. 2016), to the point that it has long been debated whether the Sun might actually have originated from NGC 2682 (see Pichardo et al. 2012, and references therein). Evidence of substantial dynamical evolution in this cluster has been reported by several authors, such as mass segregation (e.g. Fan et al. 1996; Balaguer-Núñez et al. 2007; Bonatto & Bica 2003; Davenport & Sandquist 2010; Geller et al. 2015; Gao 2018) or an elongated halo that is roughly aligned with the cluster proper motion (Davenport & Sandquist 2010). Dynamical simulations of this cluster suggest that it initially was ten times more massive (Hurley et al. 2005).

The astrometric information provided by the second data release of the Gaia mission, Gaia DR2 (Gaia Collaboration 2016, 2018b), has proven a valuable tool for investigating the extra-tidal regions of several OCs in the solar neighbourhood. Tidal tails have been revealed in the Hyades (Melotte 25 Meingast & Alves 2019; Röser et al. 2019), Coma Berenices (Melotte 111 Fürnkranz et al. 2019), and Praesepe (NGC 2632 Röser & Schilbach 2019) open clusters, all of them with ages younger than 1 Gyr (e.g. Bossini et al. 2019). However, hints of the existence of cluster members beyond the tidal radius have also been reported around the cluster Ruprecht 147 (Yeh et al. 2019), which is older than 1 Gyr (e.g. Kharchenko et al. 2013).

In this paper we study the current dynamical status of NGC 2682 using Gaia DR2. The observational material used in this paper is presented in Sect. 2. The obtained density profile is discussed in Sect. 3. The 3D spatial distribution is derived in Sect. 4 using a Bayesian approach to parallax inversion. Finally, in Sect. 5 we draw conclusions, showing that NGC 2682 is more extended than expected and that its outer envelope of extra-tidal stars may be the result of disc-shocking.

2. Data

Gaia DR2 provides a five-parameter astrometric solution, that is, positions (α and δ), proper motions (μ_α*, μ_δ), and parallaxes (ϖ) plus magnitudes in three photometric bands (G, G_BP, and G_RP) for more than 1.3 billion sources, and radial velocities (RV) for more than 7 million stars (Gaia Collaboration 2018b; Katz et al. 2019; Lindegren et al. 2018). We study the data from Gaia DR2 in a radius of 10 deg (∼150 pc) from the centre of NGC 2698 as determined by Cantat-Gaudin et al. (2018a): α_c = 132$\stackrel{°}{.}$846 and δ_c = 11$\stackrel{°}{.}$814. We followed the recommendations by Arenou et al. (2018), Lindegren et al. (2018, Eqs. (C.1) and (C.2)) and the Gaia DR2 known-issues cosmos webpage¹ for an astrometrically clean subset. We also limited our analysis to objects with low parallax uncertainties. The uncertainty in parallax depends on the signal-to-noise ratio, which is correlated with G magnitude, but also on the fact that the star is a binary. Therefore we rejected stars with errors in parallaxes larger than twice the typical error at that G magnitude. Finally, we constrained our analysis to stars with a G magnitude brighter than 18 mag to avoid problems with the completeness of the sample at fainter magnitudes. We did not attempt to use the radial velocities in our analysis because Gaia only provides radial velocities for the very bright objects (G ≤ 12, see Katz et al. 2019; Soubiran et al. 2018), and other studies available in the literature (e.g. Geller et al. 2015) cover areas smaller than we explore here.

We used two different methods to determine the membership probabilities: the UPMASK (unsupervised photometric membership assignment in stellar clusters; Krone-Martins & Moitinho 2014) and Clusterix² (Balaguer-Núñez et al. 2017) codes. Although UPMASK was originally developed to handle photometric data, it can also be applied to astrometric data (see the discussion in Cantat-Gaudin et al. 2018b). UPMASK assigns the membership probabilities on the basis that member stars must share common properties in μ_α*, μ_δ, and ϖ space and are more concentrated in α − δ plane than a random distribution (see Cantat-Gaudin et al. 2018b, for details). On the other hand, Clusterix performs an empirical determination of the field and cluster frequency functions in proper motions (μ_α*, μ_δ) by comparing two areas: the cluster core, and an area far enough to represent the field population in the cluster area (see Appendix A for details). For UPMASK, we considered as cluster members stars with an astrometric probability p ≥ 0.4 (Cantat-Gaudin et al. 2018a; Soubiran et al. 2018). In the case of Clusterix, we selected as member candidates objects with p ≥ 0.81 following the procedure described in Appendix A, applied after an additional constrain in ϖ. Clusterix provides slightly more members than UPMASK. Because there is no clear reason to prefer one method over the other, we selected the objects for which both methods agree and considered them as members. Finally, we fine-tuned the sample on the basis of the positions in the colour–magnitude diagram: main-sequence (MS), turn-off (TO), sub-giant branch (SGB), red giant branch (RGB), and red clump (RC). Basically, we used objects inside a radius of 25′ to define the regions we used for this purpose (Fig. 1). After all this procedure, our sample consists of 808 objects inside a radius of 150 pc. Following the same procedure, we selected stars redward of the main sequence in the expected location of the binary sequence. There are 144 stars in total in this region (grey points in Fig. 1).

Fig. 1. Gaia DR2 colour–magnitude diagram of the stars that are selected as members by both UPMASK and Clusterix. The different populations are plotted in different colours: RGB/RC (red), SGB (magenta), TO (green), MS (purple), BSS from Clusterix (BSSC, blue squares) and UPMASK (BSSu, cyan), and binaries (grey).

Several stars also lie above the main-sequence turn-of. They are known as blue straggler stars (BSSs, see Fig. 1). The formation of these objects is still not well understood, but the proposed scenarios include mass transfer or mergers in close binary systems (e.g. McCrea 1964) or collisions between single stars or multiple systems that led to a merger (e.g. Hills & Day 1976). They are key parts for investigating the internal dynamics because they are typically located in the central regions, probably as a consequence of dynamical friction (e.g. Mapelli et al. 2004). According to Geller et al. (2015) and references therein, this cluster contains 25 BSS candidates (see Table 1). One of these stars, S 1082, is not present in the Gaia DR2 catalogue. This is a triple system that includes an RS CVn-type eclipsing binary (Belloni et al. 1998). Both Clusterix and UPMASK methods yield very low membership probabilities for nine of these stars. Two of them have been discarded on the basis of their radial velocities, and the rest are binaries according to Geller et al. (2015). Finally we obtained a list with 11 BSS candidates with high membership probabilities for the two methods: two of them are new findings of this work in the external regions (see Sect. 3), 8 have previously been studied by Geller et al. (2015), and the last star, S 1053 is a double-lined spectroscopy binary according to Geller et al. (2015, see Table 1). Seven of the objects in common with Geller et al. (2015) were classified as “single members”, 2 were marked as “binary likely member”, one is the double-lined spectroscopy binary S 1053. The remaining 8 BSS candidates have a high membership probability for Clusterix, but not for UPMASK. Five of them are labelled by Geller et al. (2015) as “binary member” or “binary likely member”. In principle, this can be explained by the fact that UPMASK takes the uncertainties in the Gaia DR2 astrometric parameters into account, which in principle may be larger for binaries. There is no additional information for the other 3 previously not studied objects in the external regions of the cluster. Because there is no additional reason to prefer one over the other, we analyzed the Clusterix and UPMASK BSS samples independently (BSSc and BSSu, respectively).

The stellar sample used in this work, including stellar identification, coordinates, and membership probabilities derived by both UPMASK and Clusterix, is only available at the CDS.

3. Radial density profile

The radial density profile is a basic tool for investigating the spatial distribution of the cluster stellar populations and their extension. To do this, we calculated the mean stellar surface density in concentric rings as ${\rho }_i=N_i/\pi (R_{i+1}^2-R_i^2)$, where N_i is the number of stars in the i-th ring with inner and outer radius R_i and R_i + 1, respectively. For this purpose, we used only the stars along the cluster sequence, excluding BSS candidates and binaries. The density uncertainty in each ring was estimated assuming Poisson statistics.

The radial density profile is shown in Fig. 2. The stellar density decreases slowly with radius in the inner ≈30′. After this, there is a steeper decrease until ≈150′. From there on, the slope slows down but does not flatten because there are still members at larger distances, which are not shown in the figure for clarity. The first conclusion is that there are member stars at distances larger than 100′ (25 pc) from its centre. A preliminary analysis of the Gaia DR2 sample showed that the cluster extends at least beyond ∼60′ (≈15 pc) from its centre (Gaia Collaboration 2018a). The value found here is much higher than previous determinations in the literature. From 2MASS photometry, Bonatto & Bica (2003) found that the cluster extends up to a distance of about 24′ (≈6 pc), while Davenport & Sandquist (2010) used Sloan Digital Sky Survey (SDSS) photometry and reported that the cluster extends to ∼60′ (≈15 pc). More recently, Gao (2018) found a limiting radius of ∼62′(≈16 pc) also using Gaia DR2, but with a different method for assigning membership probabilities.

Fig. 2. Stellar radial density profile for M 67 members (black filled circles). The different templates are plotted in different colours (see text for details). Vertical grey lines mark the derived rc and rt for the King (short dashed lines) and Nuker (long dashed lines) templates and rc for the EFF template (dot-dashed line). Small arrows in the bottom show previous cluster size determinations by Bonatto & Bica (2003), Davenport & Sandquist (2010), and Gao (2018) from left to right.

In order to characterize the resulting density profile, we fitted it with three different analytic templates called King (King 1962), EFF (Elson et al. 1987), and Nuker (van der Marel & Anderson 2010). Although these expressions were obtained empirically and are not the physical solutions of the collisionless Boltzmann equation, they are useful for analysing the radial density profile. It is assumed that the density profile can be separated into three parts: core, bulk, and halo (or tidal debris). The different templates were defined to focus on each of these parts. For example, King (1962) tried to reproduce the bulk of the cluster, whereas Elson et al. (1987) focused on the external tidal regions.

In the following we present the analytical templates of the density profiles we investigated. A more detailed explanation of each of them is provided by Küpper et al. (2010).

• 1.

  The widely used King template was presented by King (1962). It is in the form

  $$\begin{array}{c}\hfill \rho (r)={\rho }_{\mathrm{bg}}+{\rho }_0{[{\displaystyle \frac{1}{\sqrt{1+{(r/r_{\mathrm{c}})}^2}}}-{\displaystyle \frac{1}{\sqrt{1+{(r_{\mathrm{t}}/r_{\mathrm{c}})}^2}}}]}^2,\end{array}$$(1)

  where ρ_bg is the background density and ρ₀ is the central density. The core radius r_c is defined as the distance between the centre and the point where ρ(r) = ρ₀/2. The tidal radius r_t is the position where ρ(r) = ρ_bg, the background density. The King template therefore consists of a flat core and a bulk, but has no term for the tidal debris.

• 2.

  Given that the outer parts of the density profiles of several clusters were not properly reproduced by King profiles, Elson et al. (1987), EFF proposed this template, which falls like a power law without edge radius:

  $$\begin{array}{c}\hfill \rho (r)={\rho }_0{[1+{\left({\displaystyle \frac{r}{r_{\mathrm{c}}}}\right)}^2]}^{-\eta /2},\end{array}$$(2)

  where ρ₀ and r_c have the same meaning as in the case of the King template. η is the slope of the template for radii much larger than r_c, and η is the slope of the template for radii much larger than the core radius.

• 3.

  The Nuker template was initially proposed by Lauer et al. (1995) to describe the surface brightness profiles of elliptical galaxies that showed a power-law cusp towards the centre and a logarithmic decline for radii larger than r_c. An additional logarithmic slope from r_t was added by van der Marel & Anderson (2010) in order to describer the surface brightness of ω Centauri. This template is able to fit a cluster with a non-flat core, such as a core-collapse cluster or a cluster with a massive central black hole. It has the form

  $$\begin{array}{c}\hfill \rho (r)={\rho }_{\mathit{bg}}+{\rho }_0{2}^{\frac{\eta -\gamma }{\alpha }}{\left[{\displaystyle \frac{r}{r_{\mathrm{c}}}}\right]}^{-\gamma }{[1+{\left({\displaystyle \frac{r}{r_{\mathrm{c}}}}\right)}^{\alpha }]}^{-\frac{\eta -\gamma }{\alpha }}{[1+{\left({\displaystyle \frac{r}{r_{\mathrm{r}}}}\right)}^{\delta }]}^{-\frac{ϵ-\eta }{\delta }},\end{array}$$(3)

  where ρ_bg, ρ₀, and r_c are as before. γ and η are the power-law slopes inside and outside r_c, respectively. The factor α smoothes the transition between the two slopes. r_t is an additional break radii, from there, the density falls with ϵ logarithmic slope where δ determines the smoothness of the transition between ϵ and η. van der Marel & Anderson (2010) found that the two break radii roughly correspond to the core and tidal radius of the cluster. In brief, this template consists of a flexible core, a bulk, and tidal debris.

Each analytical template was fitted to the observed profile using a non-linear fitting with the MPFIT package (Markwardt et al. 2009). The best fits in each case are plotted with different colours in Fig. 2. The King template (red line in Fig. 2), consisting of a flat core and a bulk, is only able to reproduce the innermost region, ≲3′. Up to ≈60′, this template predicts fewer stars than observed. The EFF template (Elson et al. 1987), shown as a green line in Fig. 2, is unable to reproduce the central region up to ≈10′, predicting fewer stars than are observed. Conversely, for larger radii, this template reproduces the number of observed stars reasonably well. The modified Nuker template proposed by van der Marel & Anderson (2010) is in between the previous templates. It is able to fit the data with a non-flat core, such as a core-collapse cluster or a cluster with a massive central black hole. In the case of NGC 2682, this template (blue line in Fig. 2) reproduces the observed stellar density profile almost up to a radius of ≈200′. In the outermost radii, the Nuker template predicts a continuous logarithmic fall. However, the observed profile shows a change in slope, suggesting a tidal tail or a halo.

In order to verify the statistical significance of the agreement or disagreement of the templates and the observed density profile, we computed the χ² of the fit in the way

$$\begin{array}{c}\hfill {\chi }^2=\sum \frac{{({\rho }_{i,\mathrm{t}}-{\rho }_i)}^2}{{\sigma }_{{\rho }_i}},\end{array}$$

where ρ_i, t is the density predicted by template t at the i-th ring. We obtained χ² = 6825.6, 71.5, and 2.4 for the King, EFF, and Nuker templates, respectively. The values for the best King template fit are ρ₀ = 3.8 ± 0.6 stars arcmin⁻¹, r_c = 2$\stackrel{\prime }{.}$4 ± 0$\stackrel{\prime }{.}$2, and r_t = 319′ ± 9′. In the case of EFF template, the best fit is obtained for ρ₀ = 0.85 ± 0.01 stars arcmin⁻¹, r_c = 17$\stackrel{\prime }{.}$27 ± 0$\stackrel{\prime }{.}$03, and η = 4.23 ± 0.01. Finally, the values of the best Nuker template fit are ρ₀ = 6.6 ± 0.8 × 10⁻⁴ stars arcmin⁻¹, r_c = 90′ ± 2′, η = 10.2 ± 0.1, γ = 0.57 ± 0.01, α = 1.14 ± 0.01, r_t = 144′ ± 1′, δ = 28 ± 4, and ϵ = 3.1 ± 0.1. The Nuker template reproduced the observed stellar density profile for NGC 2682 better than the others, as the χ² value confirms.

Figure 3 shows the cumulative projected radial distribution of the different NGC 2682 stellar populations normalised to the total number of the sub-population and including BSS and binaries. It is clear that the global distribution (open black circles) is mainly due to the MS stars (magenta squares), particularly in the outskirts. Turn-off objects (green squares) follow a similar trend as MS stars, but they are more concentrated. While 15% of the MS population is located inside a radius r_15% of 5$\stackrel{\prime }{.}$6, the TO is concentrated inside 2$\stackrel{\prime }{.}$7 (see Table 2). r_15% increases later in the cluster sequence. The same trend is observed for r_85%, the radius that contains 85% of the stars of a given population. The binaries (grey stars in Fig. 3) are slightly more concentrated than the MS objects, as expected because they may be more massive objects. Because the mass of the star decreases with cluster sequence, the observed trend confirms that this cluster is mass-segregated; the most massive objects are concentrated in the central regions. This behaviour has been reported in the literature by several authors (e.g. Balaguer-Núñez et al. 2007; Geller et al. 2015; Gao 2018).

Fig. 3. Cumulative projected radial distribution of different populations identified in NGC 2682 normalised to the total number of stars in each population. Errorbars have not been plotted for clarity.

The distribution of the BSS deserves a more detailed analysis. In principle, the Clusterix, blue upward triangles in Fig. 3, and UPMASK, cyan downward triangles in Fig. 3, samples both show a similar trend. The bulk of the BSS are confined in the central regions of the clusters with 50% of the population inside r_50% ∼ 5′ (see Table 2). Part of the BSSs is located at large radii, although they need to be confirmed as BSS objects. Between ≈10′ and 100′ lie are almost no BSSs. The trend that the majority of the objects lies in the central region and several lie in the outskirts and almost none lie in between has been reported in globular clusters (e.g. Ferraro et al. 1997; Lanzoni et al. 2007). The observed distribution could be due to the BSS formation mechanism, as suggested by dynamical simulations (Mapelli et al. 2004). The objects in the external regions would be the result of mass transfer in primordial binaries where BSSs in the cluster core have most likely a collisional origin.

4. 3D spatial distribution

In order to further investigate the outskirts of NGC 2682, we derived its 3D spatial distribution. Without observational uncertainties on the parallaxes, the individual distances could be simply obtained as the inverse of the parallax, d = ϖ⁻¹. However, observational errors on measured parallaxes are not negligible, and their direct propagation to distances through inversion may yield unacceptable results that would translate into a strong elongation of the 3D shape along the line of sight (Smith & Eichhorn 1996; Bailer-Jones 2015; Luri et al. 2018). To overcome this difficulty, we used a Bayesian approach to parallax inversion that we describe in Appendix B, essentially following Bailer-Jones (2015), but with a different prior that combines two terms: the cluster proper, or core, and a field or halo term that accounts for stars that either form a diffuse envelope around the cluster or are serendipitously placed at a similar distance and share proper motions with the cluster. The first term is represented by a normal distribution for simplicity, and the second term is the exponentially decreasing volume-density prior introduced by Bailer-Jones (2015). A detailed analysis of the prior that is best for clusters is beyond the scope of this paper. In our case, the form of the prior does not have a strong effect because the main result comes without 3D information, that is, the sky-projected positions of stars that does not depend on parallaxes, but shows a small tail in the opposite direction to that of the movement (top right panel of Fig. 4). On the other hand, if the parallaxes yield enough information, then the shape of prior makes no great difference unless it is very unrealistic.

Fig. 4. Position of the M 67 members in Galactic Cartesian coordinates (grey dots). Contours show different stellar densities. Red arrows are proportional to the velocities in each axis derived by Soubiran et al. (2018). Blue arrows show the direction to the observer.

The distribution of NGC 2682 in Galactic Cartesian coordinates is shown in Fig. 4. In this coordinate system, the X axis points towards the Galactic centre, Y is positive towards the local direction of rotation in the plane of the Galaxy, and Z points towards the North Galactic Pole. The Galactic centre is the origin of this system, with the Sun located at (−8340,0,27) pc (see Katz et al. 2019, and references therein). The spatial distribution of the stars is not isotropic in the XZ- and ZY-planes, although the symmetric nature of the normal term represents the cluster in the prior. This is a clear indication that the observed parallaxes have a non-negligible effect on the posterior. The XZ-plane (top left panel) shows a clear tail almost in the direction opposite to that of the cluster movement. An elongation is observed in the direction of the movement in the ZY-plane (bottom right panel). A small elongation is also observed in the XY-plane (bottom left panel), but in this case, in the direction perpendicular to the movement.

Davenport & Sandquist (2010) used ground-based photometry and reported an asymmetric halo up to a distance of 60′ from the cluster centre that is elongated roughly in line with the proper-motion vector on the sky plane (see Fig. 10 of Davenport & Sandquist 2010). However, no such elongation is observed in our case (top right panel of Fig. 4). However, as commented before, a small tail is visible in the opposite direction to that of the movement.

5. Physical interpretation

We estimated the radius of the NGC 2682 Hill sphere R_H in the Galactic potential by approximating the latter to a point-particle potential. We obtained

$$\begin{array}{c}\hfill R_{\mathrm{H}}=D\times {\left(\frac{m}{3M}\right)}^{1/3},\end{array}$$(4)

where D is the distance of NGC 2682 from the Galactic centre, m is the mass of NGC 2682, and M is the mass of the Galaxy enclosed within the orbit of NGC 2682. The quantity we know with less uncertainty is D, which is estimated as D = 8.9 kpc from Gaia DR2 (Cantat-Gaudin et al. 2018a) On the other hand, the total mass of the Galaxy is model dependent, for instance, M = 10¹² M_⊙ (McMillan 2017), and this uncertainty is compounded by the fact that we wish to consider only the mass within the orbit of NGC 2682. Similarly, there is likely an uncertainty of a factor 2 in estimates of the NGC 2682 mass such as Hurley et al. (2005), who set it at m = 2 × 10³ M_⊙. Fortunately, the dependence to the power of 1/3 in Eq. (4) ensures that the uncertainty on these numbers has little effect on R_H.

To illustrate this, we consider two different extreme scenarios:

• Overestimate the enclosed Galactic mass at M = 10¹² M_⊙ and underestimate the NGC 2682 mass at m = 10³ M_⊙ to obtain R_H = 6.2 pc

• Underestimate the enclosed Galactic mass at M = 10¹¹ M_⊙ and overestimate the NGC 2682 mass at m = 10⁴ M_⊙ to obtain R_H = 28.6 pc

In both cases, our observations show that M 67 stars can be found farther out than R_H, suggesting that they are indeed extra-tidal. Such loosely bound stars may have originated as the consequence of internal relaxation (i.e. they may be dynamically evaporated stars) or as the effect of an external disturbance.

We used the software GRAVPOT16 (Fernandez-Trincado 2017; Robin et al. 2012, 2003) to integrate the orbit of NGC 2682 back in time for 1 Gyr. We find that the last disc passage of NGC 2682 occurred ≈40 Myr ago, as shown in Fig. 5, which plots the height over the Galactic plane Z as a function of time in the last 200 Myr. The shadowed region shows the width of the thin disc according to Jurić et al. (2008). The thick-disc width is larger than the plotted area. This result is qualitatively compatible with the findings of Pichardo et al. (2012), even though the exact time of the last crossing differs due to systematics introduced by the usage of different Galactic potential models and to the lack of good proper motions in the pre-Gaia era.

Fig. 5. Height on the Galactic plane Z as a function of time in the last 200 Myr of the orbit of NGC 2682 as integrated by GRAVPOT16. The green star shows the time of the last disc passage ≈40 Myr ago. The shadowed area shows the extension of the thin disc.

According to Hurley et al. (2005), the relaxation time of a realistic N-body model of NGC 2682 was initially (before mass loss) ≈300 Myr. Even when we allow that mass loss has shortened the relaxation time by a factor ten, it is reasonable to assume that the effects of this disc passage have not yet been erased by dynamical relaxation, especially in the external regions of the clusters where the relaxation time is shorter than the half-light relaxation time. According to Fig. 5, NGC 2682 has had three disc-crossing passages in the last 200 Myr, which implies that the cluster is constantly disturbed and will never be relaxed because of the frequent passages.

The extra-tidal features we observe in M 67 would thus be quite a natural result of disc-shocking. This is to be confirmed quantitatively with dedicated N-body models in an upcoming paper (Pasquato et al., in prep.).

An estimate of the effects of disc crossing can be obtained in the analytical approximation of an impulsive shock (Ostriker et al. 1972), where every star receives a nudge in velocity,

$$\begin{array}{c}\hfill \delta v_z=\frac{4\pi G\mathrm{\Sigma }z}{V},\end{array}$$(5)

where V is the velocity with which the cluster crosses the disc, z is the distance along the Z coordinate of a given star from the cluster centre, and Σ is the density per unit surface of the Galactic disc. For a typical star that moves with a velocity v = GM/R_h relative to the cluster centre, this corresponds to a relative change in energy of order

$$\begin{array}{c}\hfill \frac{\delta E}{E}=\frac{{\delta v_z}^2}{v^2}\approx \frac{16{\pi }^{2}G{\mathrm{\Sigma }}^2{R}_{\mathrm{h}}^3}{M{V}^2},\end{array}$$(6)

where R_h is the cluster half-mass radius, and M is the total mass of the cluster. This is proportional to $M_{\text{D}}^2{\sigma }^2/M^2{V}^2$, where M_D is the mass enclosed on the Galactic disc by a circumference with a radius equal to the cluster half-mass radius, and σ is the velocity dispersion of the cluster.

Based on our orbital calculations, V ≈ 30 km s⁻¹, and a surface mass density of 30 M_⊙ pc⁻² (McGaugh 2016), we obtain δE/E ≈ 0.3, which corresponds to an expansion of the virial radius of the cluster of the same quantity, that is, δR/R ≈ δE/E ≈ 0.3. This is enough to justify a relatively large number of extra-tidal stars.

6. Conclusions

Two different methods, Clusterix and UPMASK, were used to assign membership probabilities to stars in the line of sight of NGC 2682 up to a distance of 150 pc (10°) using the information provided by the Gaia-DR2 catalogue. Additionally, we estimated distances to obtain 3D stellar positions using a Bayesian approach to parallax inversion, with an appropriate prior for star clusters. The main results of our analysis are listed below.

• The cluster extends up to 200′ (50 pc), which implies that its size is at least twice as large as previously believed. This exceeds the cluster Hill sphere based on the Galactic potential at the distance of NGC 2682.

• Only a template with a power-law cusp towards the centre and a logarithmic decline for larger radii is able to reproduce the observed stellar density profile. In any case, there are stars outside the tidal radius.

• As reported in previous studies, we confirm that the cluster is mass-segregated; the most massive objects are concentrated in the central regions.

• As observed in several globular clusters, the BBSs show a bimodal distribution in which the majority of the objects are located in the core an several objects lie in the outskirts, with nothing in between. The objects at larger radii remain to be confirmed as real BSSs, however.

• The extra-tidal stars in NGC 2682 may originate from external perturbations such as disc-shocking or dynamical evaporation from two-body relaxation. The former origin is plausible given the orbit of NGC 2682: the cluster crossed the Galactic disc ≈40 Myr ago. We plan to further investigate this scenario with direct N-body models.

Acknowledgments

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 664931 This work was partly supported by the MINECO (Spanish Ministry of Economy) through grant ESP2016-80079-C2-1-R (MINECO/FEDER, UE) and MDM-2014-0369 of ICCUB (Unidad de Excelencia “María de Maeztu”), European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement GENIUS FP7 – 606740 and the European Commission Framework Programme Horizon 2020 Research and Innovation action under grant agreement ASTERICS 653477. Support from PREMIALE 2016 MITiC is acknowledged. Based on Clusterix 2.0 service at CAB (INTA-CSIC). DB is supported in the form of work contract FCT/MCTES through national funds and by FEDER through COMPETE2020 in connection to these grants: UID/FIS/04434/2019; PTDC/FIS-AST/30389/2017 and POCI-01-0145-FEDER-030389

References

Appendix A: Clusterix membership selection

Clusterix is a web-based tool for the determination of membership probabilities from proper-motion data using a non-parametric approach based on the formalism described by Galadi-Enriquez et al. (1998), in a new implementation whose details can be found in Balaguer-Núñez et al. (2017). In the sky area occupied by the cluster, the frequency function is assumed to consist of two contributions: cluster and field stars. The tool performs an empirical determination of the frequency functions from the vector-point diagram (VPD) without relying on any previous assumption about their profiles. To distinguish the two populations, Clusterix studies the VPD corresponding to two areas on the plane of the sky: one in the cluster core, and another one far enough from the cluster to have a small (negligible) cluster contribution, but close enough to still be a good representation of the field in the cluster area. Clusterix allows searching the appropriate areas in an interactive and visual way until an optimal separation is obtained (more details will be given in Balaguer-Núñez et al., in prep.). Several limits can be adjusted to make the calculation computationally feasible without interfering with the quality of the results: in this case, we have set a constraint on proper motions (−15 ≤ μ_α* ≤ +15 and −15 ≤ μ_δ ≤ +15 mas yr⁻¹), proper-motion errors (below 5 mas yr⁻¹), and magnitude (G <  18). Because the Clusterix selection is only based on proper motions, we added an additional constraint in parallaxes by choosing only objects within ±2σ from the median of the ϖ distribution. The current version of Clusterix does not apply any criteria based on spatial distribution or photometry for membership assessment.

There is no rigorous way to decide where to set the limit between members and non-members in a list sorted by membership probabilities. This non-parametric membership method produces an expected number of cluster members that serves as an indication of the probability limit for the most probable cluster members. However, when a large sky area is studied, as we did here, the field surrounding the cluster cannot be considered perfectly homogeneous because other structures exist in the background and surroundings (e.g. other clusters). This is not a problem for the membership calculation as long as the area chosen to represent the cluster is not too far from the cluster core, but it has an effect on the calculation of the expected number of cluster members because this quantity depends on the total number of stars in the area. In this case, we have performed this study in two steps: first considering an area of 5 deg where the field is homogeneous enough and the number of expected cluster members gives a reliable result, followed by a second step, using that list of members (where most of the cluster is found) to fix the membership probability cut on the study of the 10 deg area (p ≥ 0.81).

Appendix B: Bayesian inference of stellar distances for reconstructing the 3D shape of NGC 2682

As commented in Sect. 4, the propagation of the non-negligible ϖ uncertainties causes the determination of distances through the inversion of the parallax, d = ϖ⁻¹ to produce erroneous results. In order to overcome this problem, we have used a Bayesian approach to determine distances d to the objects we analysed. To do this, we followed the Bailer-Jones (2015) Bayesian approach and used for each star the parallax ϖ to update a prior distribution in distance Π(d). We obtained an (un-normalised) posterior distribution as follows:

$$\begin{array}{c}\hfill f(d|\varpi )\propto f(\varpi |d)\mathrm{\Pi }(d),\end{array}$$(B.1)

where

$$\begin{array}{c}\hfill f(\varpi |d)\propto {\mathrm{e}}^{-{(\varpi -\frac{1}{d})}^2/2{\sigma }_{\varpi }^2},\end{array}$$(B.2)

with ${\sigma }_{\varpi }^2$ equal to the quoted error bar on each parallax measurement, that is, we assumed that uncertainties on the parallax are Gaussian, with a standard deviation equal to the quoted error bar σ_ϖ. For the prior, we used a Gaussian centred at the average distance of NGC 2682, d₀ ≈ 860 pc (Cantat-Gaudin et al. 2018a), to represent the cluster, plus a term that accounts for field stars, or more in general, stars that may be forming a diffuse envelope around NGC 2682, for which we do not wish to make overly restrictive assumptions. The standard deviation of the Gaussian term σ_d was chosen equal to the standard deviation of physical distances of member stars projected on the plane of the sky. This is equivalent to assuming in the prior that the cluster term is spherical. While the choice of a Gaussian prior may seem unjustified from the physical point of view, we chose it for simplicity because it is significantly easier to implement than a King model, for instance. The resulting prior is as follows:

$$\begin{array}{c}\hfill \mathrm{\Pi }(d)=\alpha \xi (d)+(1-\alpha )\frac{1}{\sqrt{2{\sigma }_d}}{\mathrm{e}}^{-{(d-d_0)}^2/2{\sigma }_{\mathrm{d}}^2},\end{array}$$(B.3)

where ξ(r) is a (normalised) exponentially decreasing volume-density prior (see Eq. (17) of Bailer-Jones 2015, with L = 8 kpc), and 1 − α is the UPMASK membership score. As the Bailer-Jones (2015) prior is essentially a smoothly truncated constant density prior, it is a good choice to avoid strong assumptions about the stars that by chance are loosely associated with NGC 2682. The adopted distance for each star is estimated with the mode of the posterior distribution. As pointed out in Luri et al. (2018), the whole posterior distribution should be kept, rather than summarised, as the end result of a correct Bayesian inference procedure. This applies in our case a fortiori because we are not interested in the position of individual stars as much as we are interested in the overall mass density distribution of the whole cluster. However, for the sake of concreteness, we here decided to obtain a point estimate from the posterior distribution.

All Tables

All Figures

	Fig. 1. Gaia DR2 colour–magnitude diagram of the stars that are selected as members by both UPMASK and Clusterix. The different populations are plotted in different colours: RGB/RC (red), SGB (magenta), TO (green), MS (purple), BSS from Clusterix (BSSC, blue squares) and UPMASK (BSSu, cyan), and binaries (grey).
In the text

Fig. 2. Stellar radial density profile for M 67 members (black filled circles). The different templates are plotted in different colours (see text for details). Vertical grey lines mark the derived rc and rt for the King (short dashed lines) and Nuker (long dashed lines) templates and rc for the EFF template (dot-dashed line). Small arrows in the bottom show previous cluster size determinations by Bonatto & Bica (2003), Davenport & Sandquist (2010), and Gao (2018) from left to right.

In the text

	Fig. 3. Cumulative projected radial distribution of different populations identified in NGC 2682 normalised to the total number of stars in each population. Errorbars have not been plotted for clarity.
In the text

	Fig. 4. Position of the M 67 members in Galactic Cartesian coordinates (grey dots). Contours show different stellar densities. Red arrows are proportional to the velocities in each axis derived by Soubiran et al. (2018). Blue arrows show the direction to the observer.
In the text

	Fig. 5. Height on the Galactic plane Z as a function of time in the last 200 Myr of the orbit of NGC 2682 as integrated by GRAVPOT16. The green star shows the time of the last disc passage ≈40 Myr ago. The shadowed area shows the extension of the thin disc.
In the text