The longevity of the oldest open clusters: Structural parameters of NGC 188, NGC 2420, NGC 2425, NGC 2682, NGC 6791, NGC 6819

Context. Open clusters’ dynamical evolution is driven by stellar evolution, internal dynamics and external forces, which according to dynamical simulations, will evaporate them in a timescale of about 1 Ga. However, about 10% of the known open clusters are older. They are special systems whose detailed properties are related to their dynamical evolution and the balance between mechanisms of cluster formation and dissolution. Aims. We investigate the spatial distribution and structural parameters of six open clusters older than 1 Ga in order to constrain their dynamical evolution, and longevity. Methods. We identify members using Gaia EDR3 data up to a distance of 150 pc from each cluster’s center. We investigate the spatial distribution of stars inside each cluster to understand their degree of mass segregation. Finally, in order to interpret the obtained radial density proﬁles we reproduced them using the lowered isothermal model explorer with PYthon ( LIMEPY ) and spherical potential escapers stitched ( SPES ). Results. All the studied clusters seem more extended than previously reported in the literature. The spatial distributions of three of them show some structures aligned with their orbits. They may be related to the existence of extra tidal stars. In fact, we ﬁnd that about 20% of their members have enough energy to leave the systems or are already unbound. Together with their initial masses, their distances to the Galactic plane may play signiﬁcant roles in their survival. We found clear evidences that the most dynamically evolved clusters do not ﬁll their Roche volumes, appearing more concentrated than the others. Finally, we ﬁnd a cusp-core dichotomy in the central regions of the studied clusters, which shows some similarities to the one observed among globular clusters.

ing gas-free clusters is driven by relaxation.The stars randomly exchange energy via gravitational interactions, causing equipartition of energy between stars of different masses.This causes a mass segregated system with the most massive objects concentrated in the centre while the lower mass stars migrate to the outskirts, forming a halo.Some of these stars acquire enough velocity to escape from the system, resulting in its gradual evaporation (e.g.Pang et al. 2021).This dissolution is amplified by the forces acting on these systems as they orbit in the Galaxy, such as encounters with giant molecular clouds or passes through the disc (e.g.Binney & Tremaine 2008;Gieles et al. 2006).Therefore, according to dynamical simulations a typical cluster, with 10 4 M ⊙ , will evaporate in a timescale of ∼1 Ga (e.g.Baumgardt & Makino 2003;Lamers et al. 2005a;Bastian & Gieles 2008).
Between 8 and 10% of the more than 6,800 known and candidate open clusters have ages older than 1 Ga (Cantat-Gaudin et al. 2020;Hunt & Reffert 2023).The properties of this older population are related to the dynamical evolution of clusters and the balance between mechanisms of cluster formation and dissolution (Friel 2013).These clusters are preferentially found at galactocentric distances larger than 6 kpc, and at larger heights from the Galactic plane.The majority are found close to their maximum excursion from the plane, where they spend most of the time, away from the disc disruptive influence.
The oldest systems are typically larger than intermediate-age clusters, 50 Ma to 1 Ga (e.g.Friel 2013).This, together with the preferential location in the outer disc, could be interpreted as larger clusters easily survive at larger distances.Nevertheless, this could be due to a selection effect, since small clusters are more difficult to detect.In spite of this, the longevity of the oldest open clusters is still not well understood.
The goal of this paper is to determine structural parameters and study the spatial distribution of stellar populations inside six of the oldest open clusters: NGC 188, NGC 2420, NGC 2425, NGC 2682, NGC 6791, and NGC 6819.This paper is organised as follows.The observational material used and astrometric membership probability determination are described in Sect. 2. The radial density profiles and their interpretation based on dynamical models are presented in Sect.3. The segregation in mass of studied clusters is investigated in Sect. 4. Several physical parameters such as Jacobi radius, half-mass relaxation time or initial mass are estimated in Sect. 5.The results are discussed in the context of the open cluster dynamical evolution in Sect.6.Finally, the main conclusions of this study are listed in Sect.7.

Membership determination
Our analysis is based on the positions (α, δ), proper motions (µ α * , µ δ ) and parallaxes (̟) provided by Gaia EDR3 (Gaia Collaboration et al. 2021;Lindegren et al. 2021) alongside magnitudes in three photometric bands G, G BP and G RP (Riello et al. 2021).We applied the G-band corrections recommended by the Gaia team 1 .Note that these are the most updated values for astrometric and photometric magnitudes provided by Gaia since DR3 only propagated the values already released in EDR3 for them (Gaia Collaboration et al. 2022).We performed our analysis in a radius of 150 pc from each cluster centre, as determined by Cantat-Gaudin et al. (2020) and listed in Table 1.This value was selected because it was larger than the known cluster sizes at that time.In order to reduce the size of the sample, we applied weak constraints in proper motions and parallaxes, discarding stars outside five times the uncertainties for proper motions and parallaxes, and centred on their average values, again using the results reported by Cantat-Gaudin et al.
(2020, see Table 1).Moreover, we limited our sample to stars brighter than G=18.5 mag to ensure a good completeness of the sample, with reasonable average uncertainties in proper motions and parallaxes.The central regions of each cluster were used to constrain their sequences in the colour-magnitude diagram, removing those objects which were far away from the cluster sequence or the position of the blue straggler stars.
For each star in these initial samples, the probability of belonging to each cluster was determined from its proper motions and parallaxes using UPMASK (Unsupervised Photometric Membership Assignment in Stellar clusters, Krone-Martins & Moitinho 2014).This tool, originally developed to assign membership probabilities from the photometric data, was adapted to do so by using astrometric data (Cantat-Gaudin et al. 2018b;Pera et al. 2021).UPMASK uses a k-means clustering algorithm assuming that the member stars are closely clustered together in µ α * , µ δ , and ̟ 3D space (see Cantat-Gaudin et al. 2018b, for a full discussion of this).Initially, we planned to follow the same procedure used by Carrera et al. (2019b) for NGC 2682 (M67).This procedure works well at large radii for clusters with average proper motion and parallax values statistically different from the average values of field stars, such as NGC 2682.However, in the cases where this does not happen, UPMASK works well at short radii where the clusters' stars dominates.If a large radius is used, UPMASK assigns lower membership probabilities even to the central members in comparison with the values obtained using a short radius.To overcome this issue, the initial sample containing stars within the 150 pc radius was split into multiple data sets containing the objects within increasing radius values.The initial radius was selected as the one which contains 50% of the cluster members, as determined by Cantat-Gaudin et al. (2020).If the sample has a projection on sky larger than 3 • .0 the data sets would have radii values increasing in steps of 0 • .5.Alternatively, if the radius of the sample is less than 3 • .0 the radius for each data set would increase by 0 • .2 each time.UPMASK was run in each data set independently.This means that the objects in the central region of each cluster have multiple membership probabilities determinations, while the stars in the outermost radius have only one.For those objects with multiple determinations, we simply assumed as membership probability the maximum value obtained, which typically is derived in the innermost radii in which this object is sampled.This provided a more exhaustive view of the centre of the cluster, whilst also keeping the maximum amount of members towards the outskirts.Finally, we considered as cluster members those objects with membership probability, p, greater or equal than 0.4 in this final merged catalogue (see Soubiran et al. 2018;Carrera et al. 2019a, for details).The impact of this selection in our results are discussed in Sect. A.
In the case of NGC 2682 we found a total of 1 170 objects with p ≥0.4 using the Gaia EDR3 data.This number is slightly higher than the number of stars found by Carrera et al. (2019b) of 1 149 from Gaia DR2 data.Not all the objects in the Carrera et al. (2019b) sample were recovered, and instead other stars appeared with high membership probabilities.This is explained by the change in the proper motions, parallaxes, and above all, in their related uncertainties between Gaia DR2 used by Carrera et al. (2019b) and EDR3 used here.Moreover, there is no preferential spatial location for both the new recovered and the discarded stars.Also, this issue only affects a small fraction of objects, ∼5%.Therefore, this ensures that the procedures adopted here and by Carrera et al. (2019b) are equivalent.
In the case of NGC 6819, Cantat-Gaudin et al. (2018a) reported 1 715 objects in a radius of 0 • .45 with p ≥0.4 and G≤18.5 mag.We find only ∼1 300 stars brighter than G=18.5 mag and with p ≥0.4 if we run UPMASK in the whole 3 • .3radius (150 pc at the distance of this cluster).With the procedure adopted here of performing the analysis in different increasing radii, we recover 1 785 objects with G≤18.5 mag and p ≥0.4.All these numbers were obtained before applying constraints in the positions of the stars in the colour-magnitude diagram.
The limited capabilities of Gaia for observing dense areas, such as the centre of the clusters, can affect the completeness of our analysis.According to Fabricius et al. (2021) the completeness of the sample is reduced for stars fainter than G∼19 mag for regions with densities around 10 5 stars deg −2 .Taking into account the distance of the clusters in our sample, only the faintest objects in their very central regions may be affected by this effect, except for NGC 2682, which is close enough to avoid problems with the crowding.
Gaia DR3 provides radial velocities for stars with G<14 mag (Gaia Collaboration et al. 2022;Katz et al. 2022).We used this information to evaluate the contamination, at least within the brightest stars, by objects with discrepant radial velocities.Owing to the large uncertainties of the Gaia DR3 radial velocities at the faint end, we only consider as real non cluster members those stars which very discrepant radial velocities, larger than 15 km s −1 , from the average value provided by Tarricq et al. (2021) for each cluster.With this procedure, we remove stars only on four clusters.We removed two objects from a total of 26 stars with high astrometric memberships, in the case of NGC 2425.Five objects were rejected in NGC 2682 from a sample of 487 objects.This filter has a higher impact in NGC 2420 and NGC 6819 where 10 and 11 objects are discarded from the radial velocity sample of 69 and 124 stars, respectively.
The total number of objects for each cluster after taking into account the position of the stars in the colour magnitude diagrams and their radial velocities are listed in Table 1.For each cluster, the selected members are manually separated into different groups as a function of their evolutionary stage from the position in the colour-magnitude diagrams.These groups are mainsequence (MS), turn-off (TO), sub-giant branch (SGB), red giant branch (RGB) which includes also the red clump (RC), candidate blue straggler stars (BSS) and candidate binaries.They are plotted with different colours in Fig. 1.We include blue straggler stars and binaries in spite of the well known limitations of UPMASK to assign high membership probabilities to these objects (see Carrera et al. 2019b, for a detailed discussion).

Radial density profiles
In order to investigate the spatial distributions of the stellar populations inside the studied clusters, it is necessary to remove projection effects, which are important in cases such as NGC 188 due to its location near the North celestial pole.For this purpose, we computed projected Cartesian coordinates, x and y, using the Eq. 2 by Gaia Collaboration et al. (2018b).In our case, we selected as origin, (α 0 , δ 0 ), the cluster centre listed in Table 1 taken from Cantat-Gaudin et al. (2020).In this system, the x-axis is antiparallel to the right ascension axis, and the y-axis parallel to the declination axis.We use the distances to each cluster reported by Cantat-Gaudin et al. (2020) and listed in Table 1 to convert x and y coordinates in parsecs.Finally, the radial distance for each star was obtained from the Cartesian coordinates defined above as r = x 2 + y 2 .
The spatial distributions of the studied clusters projected in the sky in Cartesian x and y coordinates are shown in Fig. 2, where the size of the points is related to their membership probabilities.The first noticeable feature is the tail towards the southeast of NGC 188, almost on the direction of motion marked by the arrow.Although without a clear tail, NGC 2420 and NGC 2425 show non-isotropic distributions in their outskirts.In the case of NGC 2425 it seems that it is aligned with the orbit.The other three clusters have a more regular distribution, but in the case of NGC 6819 the central region seems to be elongated in the SE to NW direction, again almost aligned with the movement, as already reported by Kamann et al. (2019) The radial density profile is the basic tool used to investigate the spatial distribution of stellar systems, such as open clusters, providing clues about their dynamical evolution.For the studied systems, we calculated the mean stellar surface density of objects within concentric rings as , where N i is the number of stars within the i-th ring with an inner radius of R i and an outer radius of R i+1 .The density uncertainty in each ring was estimated assuming Poisson statistics.
The obtained radial density profiles are shown in Fig. 3.Only two of the studied clusters, NGC 188 and NGC 2682, show some hints of flattening in the outskirts.In contrast, for the other four systems, with a change of slope, the radial density profiles continue to fall until 150 pc, which may imply that we did not reach the end of the cluster.
There are also some differences in the core regions.While some clusters show clear flat cores, such as NGC 6791, other show a cusp profile, such as NGC 2682 or in a lesser degree NGC 2420.In the case of NGC 2425, it seems that the core is smaller than in the other systems and its central region has been sampled with only two rings.The dichotomy of flat and cusp cores is well known in Galactic globular clusters (e.g.Djorgovski & King 1986;Trager et al. 1995).We will discuss this issue in detail in Sect.6.

LIMEPY models
The most direct way of interpreting the observed radial density profiles is by comparison with the prediction of dynamical models.It is widely stated in the literature that in spite of the irregular appearance of open clusters, their radial density profiles are well reproduced by isothermal and spherical King models (King 1966).The exception are the outermost external regions, which require an additional power-law decrease term (e.g.Carrera et al. 2019b).Davoust (1977) showed that the widely used King and Wilson (in the non-rotating and isotropic limit, Wilson 1975) models are particular cases of a more general family of models.These were extended to a more general class of (isotropic) lowered isothermal models by Gomez-Leyton & Velazquez (2014).Gieles & Zocchi (2015) further expanded them by introducing parametrised prescriptions for the energy truncation, related to the edge of the cluster, and for the amount of radially biased pressure anisotropy, which determines the size of the isotropic cores.They introduced the Lowered Isothermal Model Explorer in PYthon (LIMEPY). 2These models are particularly suited to describe the phase-space density of stars in tidally limited, masssegregated stellar clusters in all stages of their life-cycle.
To identify one model within the LIMEPY family, it is necessary to specify the values of five parameters: the central dimensionless potential, W 0 ; the anisotropy radius, r a ; the truncation parameter, g; the total mass of the system, M cl ; and the halfmass radius, r h .The central dimensionless potential, W 0 , is used as a boundary condition to solve the Poisson equation, and it determines the shape of the radial profiles of some relevant quantities.The anisotropy radius is related to the amount of anisotropy present in the system: the smaller it is, the more anisotropic is the model.The truncation parameter sets the sharpness of the truncation in energy: the larger it is, the more extended the models are, and the less abrupt the truncation is.When considering the isotropic version of the models, g = 0, corresponds to the Woolley (1954) models, g = 1 to the King (1966) models, and g = 2 to the non-rotating Wilson (1975) models.We added a sixth parameter, as suggested by Zocchi et al. (2016), to account for the background density, ρ bg .

SPES models
As mentioned above, LIMEPY models are not able to reproduce the outermost radii sampled in each cluster.The LIMEPY models provide a more elaborate description of stars near the truncation, but do not include the effect of the Galactic tidal potential, unlike other models (e.g.Varri & Bertin 2009).The tidal field makes the potential in which the stars move anisotropic, and it slows down the escape of stars (Fukushige & Heggie 2000;Baumgardt 2001), because escape is limited to narrow apertures around the Lagrangian points.The result is the existence of a stellar population, known as potential escapers, which is energetically unbound, but have not yet escaped because they have not reached the Lagrangian points (e.g.Daniel et al. 2017).These objects are the responsible for an elevation of the density and velocity dispersion near the Jacobi radius (Küpper et al. 2010;Claydon et al. 2017).The presence of potential escapers in globular clusters have been suggested as the responsible for the peculiarities observed in their outskirts.
These spherical potential escapers stitched models (hereafter SPES models, Claydon et al. 2019) have an energy truncation similar to the LIMEPY models discussed above, with the fundamental difference that the density of stars at the truncation energy can be non-zero.More importantly, the models include stars above the escape energy, with an isothermal distribution function that continuously and smoothly connects to the bound stars.
Apart from W 0 , M cl , and r h , the SPES models depend on two additional parameters, B and η.The value of B can be 0≤ B ≤1, where B =1 implies that there are no potential escapers.The parameter η is the ratio of the velocity dispersion of the potential escapers over the velocity scale, and it can have values 0≤ η ≤1.For η =0 there are no potential escapers.For fixed B, the fraction of potential escapers correlates with η.For a fixed η, the fraction of potential escapers anticorrelates with B, when B is close to 1.For small values of B, the fraction of potential escapers is approximately constant or correlates slightly with B at η constant.Finally, in the presence of potential escapers the SPES models are not continuous at the tidal radius, r t , but they have the ability to be solved (continuously and smoothly) beyond r t to mimic the effect of escaping stars (see Claydon et al. 2019, for a detailed discussion).

Model fits
In order to fit the observed radial density profiles with both LIMEPY and SPES models, we use a Bayesian approach following Zocchi et al. (2016) to determine the posterior probability distribution of the input parameters.For LIMEPY we choose uniform priors over the following ranges: 0.8< W 0 <15, 0.2< g <2.1, 0.1< M cl <10 6 M ⊙ ,0.2< r h <30 pc, −1< log r a <3, and −8 < log ρ bg <-1.In the case of the SPES models together with the values of W 0 , M cl , and r h values described above, we select: 0≤ B ≤1 and 0≤ η ≤1.We use a Markov chain Monte Carlo fitting technique to explore the parameter space and to efficiently sample the posterior probability distribution for the parameters above from the LMFIT PYthon implementation 3 .
We consider the widely used King, g=1, and the isotropic and non-rotating Wilson models, g=2.They provide a fairly simple description of cluster morphology, with their shape entirely determined by the dimensionless central potential W 0 .In this case, high values of W 0 implies more concentrated models.The third case is the isotropic, single-mass LIMEPY models which fit simultaneously W 0 and the truncation parameter g.Finally, we fit the SPES models.In each case, we performed 500 realisations for each cluster.

Results
The best fits for King (green dashed lines), Wilson (purple dotdashed lines), LIMEPY (blue solid lines), and SPES (red solid lines) models are shown in Fig. 3.The corresponding parameters are listed in Tables 2, 3, and 4 for King, Wilson, LIMEPY, and SPES models, respectively.Individual fits, together with their un- As expected, LIMEPY models, not only in the specific King and Wilson prescriptions but also the general ones, are not able to reproduce the outskirts of the observed radial density profiles (Fig. 3).The SPES models reproduce the profiles at large radii, due to the inclusion of potential escapers, being the ones with the lowest χ 2 values.The only exception is NGC 188 for which the fit of the LIMEPY model produced a slightly lower χ 2 than SPES ones.In any case, none of the used models are able to reproduce the cusp core observed in NGC 2682 and in a less degree in NGC 2420.On the contrary, they reproduce quite well the flat core of NGC 6791 and the intermediate region between the core and the tidal radius for all the studied systems (Fig. 3).
In general, the simple King and Wilson models produce larger χ 2 values than the general LIMEPY models without significant differences among them except in the outskirts.The smooth variation between the King and Wilson models, and also the Woolley ones, is controlled by the variation of the truncation parameter.In our case, we found that the majority of the studied clusters are close to the King models with g ∼1, contrary to what was reported by de Boer et al. (2019) for globular clusters, with values close to the Wilson template.The individual King models reproduce better the cluster profiles in four cases, while the Wilson ones work better for NGC 2425 and NGC 6819.This may be due to the fact that these systems are relatively more ex- [pc] NGC_188 4.9±0.27.3±0.36.6±0.2 606±358 28 4.2±0.47.4±0.36.6±0.2 467±339 31 NGC_2420 7.0±0.36.1±0.5 9.9±1.1 554±346 29 8.1±0.2 6.3±0.6 35.0±5.7 568±318 19 NGC_2425 8.2±0.2 4.6±0.5 17.8±2.2204±334 22 8.1±0.2 6.3±0.6 35.0±5.7 568±318 19 NGC_2682 6.3±0.2 10.5±0.5 9.3±0.5 526±337 37 6.0±0.5 10.7±2.8 9.4±1.6 339±349 40 NGC_6791 5.0±0.5 30.7±2.0 12.7±1.21.9±0.447 5.7±0.5 26.8±1.410.6±0.8 7.9±193 58 NGC_6819 6.8±0.2 22.0±1.411.5±1.053.8±280 63 6.8±0.2 25.3±1.813.9±1.7 216±337 51 [pc] [pc] [pc] NGC_188 4.7±0.40.10±0.150.62±0.026.8±0.4 6.2±0.3 3.2 17.2 0.15 0.73 37 NGC_2420 6.4±0.4 0.0±0.020.68±0.025.1±0.47.7±0.7 2.4 25.0 0.18 1.0 23 NGC_2425 7.0±0.5 0.02±0.090.74±0.033.0±0.5 9.6±2.02.1 30.5 0.23 1.2 15 NGC_2682 6.1±0.3 0.0±0.020.65±0.029.4±0.5 8.2±0.5 2.9 26.5 0.15 1.0 36 NGC_6791 4.7±0.40.06±0.20 0.65±0.0220.5±0.8 7.9±0.34.1 21.1 0.17 A key parameter in the LIMEPY models is the anisotropy parameter, r a .We found large values of r a for all the studied clusters, which implies that the amount of anisotropy in them is quite low.The only exception is NGC 6791 for which we found a very low value of r a .About the half mass radius, r h , which encloses half of the total mass of the system, the SPES models values are slightly lower than those of the LIMEPY ones with the only exception of NGC 188, for which the best LIMEPY model has a slightly lower value than in the best SPES one.In both cases, the lowest value is obtained for NGC 188.About the masses, we also found that the LIMEPY models produce slightly larger values than the SPES ones.This is explained that the SPES models consider that a fraction of observed stars are unbound to the systems, so the amount of mass needed is lower.The obtained masses may be considered as a lower limit of the real value.Owing to the fact that we have a magnitude limit, G ≤18.5 mag, we do not sample the faintest objects, which affect particularly the furthest object, NGC 6791.Using different magnitude threshold, we check that the derived masses are comparable while there are objects below the turn-off, such as one magnitude, but they decrease when the limit is around the turn-off or above it.Moreover, we excluded binaries and blue stragglers stars from this analysis.
The core radius, r c , and truncation or tidal, r t radius are not obtained from the fits, but they are computed as a function of the best results for both LIMEPY and SPES models, respectively.Derived core radii for both families of models are similar, showing the same trend that the LIMEPY models values are slightly larger than the SPES models ones, again except for NGC 188.NGC 6791 shows the largest r c , while NGC 2425 has the smallest one.Tidal radii show a similar tendency.However, the values derived for LIMEPY models are very different among clusters, from the 28.8 pc of NGC 188 to the 416 pc of NGC 6819.Only NGC 188 and NGC 2682 show a flattening in the outermost radii studied, most probably due to the algorithm not being able to properly place the end of the cluster.On contrary, the SPES models provided smaller tidal radii, between 25 pc for NGC 2420 and 55 pc for NGC 188.This implies that the objects located outside these radii are extra-tidal stars which are probably escaping from the clusters.In fact, according to SPES models between 13%, for NGC 188, and 23%, for NGC 2425, of the observed stars in each system may have enough energy to escape from them.
Other works have studied the radial density profiles of the clusters in our sample, mainly by fitting them with the analytical King profile model.The majority of them have been performed in recent years, taken advantage of the different Gaia data releases (Gao 2018;Gao & Fang 2022;Zhong et al. 2022;Angelo et al. 2023;Cordoni et al. 2023).Together with the variety of algorithms used to calculate the membership probabilities, our main difference is the studied area around the cluster centre, significantly larger than in the majority of the other cases.In general, the values for the core radius reported by these works are of the order, or slightly smaller than, the values found here.There are larger discrepancies in the tidal radius determination.This is explained in part by the small area, relative to us, covered by several of these works, and by the presence of unbound stars in the same way as the differences found above between LIMEPY and SPES models.To our knowledge, the recent work by Angelo et al. (2023) is the only one that derived masses for the studied clusters, four in common with us: NGC 188, NGC 2682, NGC 6791, and NGC 6919.Using a different approach, the masses derived for these clusters are in good agreement, within the uncertainties, with the ones obtained here from the LIMEPY models.

Mass segregation
It has been widely reported in the literature that old open clusters are segregated in mass (e.g., Mathieu 1984), with massive objects concentrated in the central region while low mass stars are dispersed in the outskirts.In spite of some evidence of primordial mass segregation in very young massive clusters (e.g.Kim et al. 2006;Stolte et al. 2006), in old open clusters this should be a direct consequence of their internal dynamics.
In order to investigate the mass segregation of the studied clusters, we have obtained the cumulative projected radial distribution of the different stellar populations normalised to the total number of the objects in each of them.In spite of the uncertainties in our membership determinations for blue stragglers and binaries, we included both populations in our analysis (see Carrera et al. 2019b, for details).According to Fig. 4, all the studied clusters are segregated in mass, with the red giant branch more concentrated than the less massive objects in the main sequence.The only exception is the central part of NGC 6791 since within the innermost 2 pc there is no clear separation among the different stellar populations.Mass segregation has already been reported in the literature for the majority of the studied clusters, mostly based on membership determination from pre-Gaia era: NGC 188 (e.g.Bonatto et al. 2005;Geller et al. 2008), NGC 2420 (e.g.Leonard 1988;Peikov et al. 2002;Paparo 1982), NGC 2682 (e.g.Balaguer-Núñez et al. 2007;Geller et al. 2015;Gao 2018), NGC 6791 (e.g.Gao 2020; Platais et al. 2011), and NGC 6819 (e.g.Kalirai et al. 2001;Kang & Ann 2002;Yang et al. 2013).In general, the global distributions are mainly due to the main sequence stars, whereas turn-off and subgiant branch objects are more concentrated.To quantify this statement, we have computed the radii, which contain the 15% and 85% of each population (Table 5).The concentration of the different populations changes from one cluster to the other.For example, in the case of NGC 2682, the more massive red giants and subgiants stars are notoriously more concentrated towards the innermost parts than the other populations.It seems that in the case of NGC 188, the turn-off objects are more concentrated than red giant and subgiants objects.This could be explained by the fact that for NGC 188 the mass difference between red giant and main sequence turn-off objects should be lower than in the case of NGC 2682 since it is about 4 Ga older.For NGC 2420, there is a non-negligible quantity of binaries in the outermost radii.Finally, the subgiant population in NGC 2425 exhibits an unexpected behaviour, flattening between ∼1 and 10 pc, but it could be and effect due to being based on a few objects.
Particularly interesting is the case of the blue stragglers, which are typically more concentrated than even the red giants.The mechanisms for the formation of blue stragglers are still not fully understood.These stars are thought to derive from normal main sequence stars that have increased in mass above a single star mass typical of the turn-off through mass transfer, mergers (e.g.McCrea 1964;Paczyński 1971), or collisions in binary systems (e.g. Hills & Day 1976).As reported by Carrera et al. (2019b), these objects show a bimodal distribution in the majority of the studied clusters, as observed in globular clusters (e.g.Ferraro et al. 1997;Lanzoni et al. 2007) and predicted by dynamical simulations (e.g.Mapelli et al. 2004).
In order to quantify and compare the degree of mass segregation among the studied clusters, we used the method proposed by Allison et al. (2009).It consists of comparing the lengths of the minimum spanning tree (MST) of the most massive stars of a cluster and a set of the same number of randomly chosen objects between all populations.An MST of a set of points is the path connecting all the points, with the shortest possible path length and without any closed loops.In a given set of points, only one MST can be drawn.In our case, we defined the mass segregation ratio (MSR) as: where l random and l massive are the average lengths of the MST of N randomly selected stars between all and massive samples, respectively.The average lengths l random,massive are calculated over 100 iterations, where at each iteration, we drew a different subsample of random stars allowing us to simultaneously calculate σ random,massive , the standard deviations of the lengths of the MST.We computed the uncertainty of Λ MSR as the square root of the quadratic sum of σ random and σ massive .
We computed the MST by using the csgraph routine implemented in the scipy PYthon module (Virtanen et al. 2020).We consider the red giant sample as representative of the massive stellar population in each cluster.After several trials, we assume N = 15 stars.In principle, a Λ MSR greater than one means that the massive population is more concentrated than a random sample, and therefore that the cluster shows signs of mass segregation.The obtained values are listed in the last column of Table 5.In all the studied clusters, the red giants are more concentrated than random selected stars.Noticeable are the larger values found for NGC 2420 and NGC 2425, which are the less massive clusters in our sample (see Table 4) and also among the youngest, where a larger mass difference between giants and unevolved stars is expected.On the other hand, the most massive systems have lower Λ MSR values.In fact, the lowest value is found for the most massive system in our sample, NGC 6791, according to the previous section results, although this is also the furthest cluster which implies that our sample only includes members in the upper MS, and therefore, we are not sampling as many low-mass stars as in other clusters.

Half-mass relaxation time
The dynamical evolution inside the clusters is related to the encounters among stars.After losing a certain amount of members, the relaxation time is the moment in which the stellar system reaches the equilibrium.In other words, this is also defined as the time required for a star to lose the memory of its initial conditions.The relaxation time is locally defined, and it can vary by several orders of magnitude in different regions of a single cluster.For this reason, the relaxation time is widely calculated in reference to the system's half-mass radius, the so-called halfmass relaxation time, t rh .
Recently, Angelo et al. (2023) have determined half-mass relaxation time for a sample which includes four systems in common with us.Except for NGC 6791, our values are larger than those derived by Angelo et al. (2023) for the four clusters in common, mainly motivated by the fact that we are using larger input clusters' masses and half-mass radii since our study covers a wider area and contains a larger number of members.

Jacobi radius
Stellar clusters do not live in isolation, but they are orbiting the Galaxy inside its tidal field.The Jacoby radius, R J , also called Hill or Roche radius, of a stellar cluster is the maximum distance inside which a star is still bound gravitationally to the system taking into account the external Galactic tidal potential.Therefore, it can be considered as a prediction of the tidal radius.According to King (1962), R J can be derived as: where R p is the perigalactic radius defined as R p = R apo (1 − e) with R apo and e being the apocentre radius and the eccentricity of the cluster orbit, and M cl and M g the masses of the cluster and the Galaxy enclosed within the orbit of each cluster, respectively.
We use the orbits parameters determined by Carrera et al. ( 2022) except for NGC 2425 for which we assume the values determined by Tarricq et al. (2021).In both cases, the clusters' orbits were computed using the galpy package (Bovy 2015).In order to determine the mass enclosed, M g , within the orbit of each system, which is model dependent, we also take advantage of the galpy package assuming the MW2014 potential for the Milky Way (see Bovy 2015, for details).
The derived Jacobi radius are listed in Table 6.Using other potentials available in the literature, such as that of McMillan (2017), yields similar values.Increasing the mass of the clusters, for example using the values derived from LIMEPY models, yields slightly larger Jacobi radius.On the other hand, an increase of the Galaxy mass enclosed inside the cluster orbit tends to reduce R J .
Alternatively, the Jacobi radius can be derived from the tidal tensor of the total potential (Renaud et al. 2011).In general, their values are lower than those derived here, except for NGC 6791.This discrepancy may be due to the fact that our masses, from SPES models, are larger than those derived by Angelo et al. (2023).In any case, neither our R J estimations nor the ones by Angelo et al. (2023) are as large as the tidal radius determined by both SPES and LIMEPY models.

Initial mass
A stellar cluster losses mass due mainly to stellar evolution and tidal disruption.Lamers et al. (2005a) showed that the fraction of the cluster initial mass, M ini , lost due to stellar evolution, (∆M) ev , in the GALEV (Galaxy Evolutionary Synthesis, Kotulla et al. 2009) models can be written as (∆M) ev M ini ≡ q ev , which can be approximated by: log q ev (t) = log age − a ev (b ev ) + c ev for t > 12.5 Ma where a ev , b ev , and c ev are coefficients which slightly depend on metallicity according to Table 1 by Lamers et al. (2005a).Metallicities are obtained from iron abundances using the values provided by Carbajo-Hijarrubia et al. in prep.except for NGC 2425 which was obtained from Randich et al. (2022).More generally, the evolution with time of the mass of a cluster which has survived the infant mortality, age ≥10 7 a, can be described as: where the first term describes mass loss due to stellar evolution and the second by disruption.
The mass decrease of a cluster can be approximated very accurately by: where masses are expressed in M ⊙ , µ ev (t) = 1 − q ev (t), and γ=0.62 according to Lamers et al. (2005a).The constant t 0 depends on the galaxy potential in which the cluster is moving on, and on the eccentricity, e, of its orbit.From Lamers et al. (2005b), it can be derived from: amb where C env,0 =810 Ma for the Milky Way (see Lamers et al. 2005b, for details).ρ amb is the ambient evaluated at the position of the cluster, expressed in M ⊙ pc −3 .It was determined for each cluster with the galpy package, assuming the MW2014 potential (Bovy 2016).
The initial cluster mass can be easily derived by manipulating Eq. 1 as: The obtained M ini values are listed in Table 6.In comparison with Angelo et al. (2023), who determined initial masses following a similar procedure, our estimations are slightly larger, except for NGC 6791, the most massive system in our sample, which is significantly smaller.Together with the different input clusters' masses, the only difference between both studies is the treatment of the Galaxy tidal field.
Finally, Baumgardt & Makino (2003) found that the disruption time of a stellar cluster can be expressed as a function of M ini as t dis = t 0 M γ ini with t dis in years.Derived values are also listed in Table 6.Obtained values are discussed in the following section.

Dynamical evolutionary stages
We discuss the dynamical stages of the studied clusters as a function of the derived structural parameters (r h , r t , r c ), relaxation times (τ rh ), Jacobi radii (R J ), location in the disc (R gc , Z), and evolution-related parameters (ages, stellar masses).Our results are summarised in the different panels of Fig. 5.In order to compare to the global trends described by open clusters, we have also plotted the results obtained by Angelo et al. (2023).These authors studied four of the clusters in our sample: NGC 188, NGC 2682, NGC 6791, and NGC 6819.We recovered slightly large structural parameters, mainly due to the fact that we are covering a large area around the clusters.In either case, the results obtained for these systems, mainly the ratios among the different radii, are compatible within the uncertainties.
The two oldest clusters in our sample, NGC 188 and NGC 6791, are also the most dynamically evolved systems according to their dynamical ages, τ rh >0.7 dex.These systems have galactocentric distances between 8 and 9 kpc from the centre of the Galaxy.They are located more than 500 pc above the [dex] [pc] [10 3 M ⊙ ] [Ga] [Ma] NGC_188 0.7±0.1 1.01±0.0832.1±1.1 26.1±2.213.4±0.844±1 NGC_2420 0.7±0.1 0.42±0.1032.4±1.39.7±0.7 7.9±0.5 30±2 NGC_2425 0-9±0.20.35±0.1027.7±2.49.5±1.84.0±0.6 54±3 NGC_2682 1.0±0.1 0.55±0.0835.0±0.9 23.8±2.1 9.5±1.038±2 NGC_6791 1.4±0.20.78±0.0841.0±1.9 70.4±4.516.3±0.826±2 NGC_6819 1.5±0.20.11±0.0740.8±1.5 34.8±1.99.6±1.412±1 Galactic plane, and therefore, less affected by the destructive influence of the disc.This is in agreement with the idea that cluster of similar ages are preferentially located outside the dense disc, which may plays a role in their survival (e.g.Friel 2013;Cantat-Gaudin et al. 2020).These two systems also have the lower r t /r c , which can be explained by the fact that they have reduced their sizes throughout their lives by losing a significant part of their outskirts.The remnants that we observe nowadays would be mainly their cores.This hypothesis is supported by the fact that they have lost more than 70% of their birth masses.Moreover, these two systems show the lower mass segregation ratio, Λ MSR .This can be a consequence of the fact that they have lost a significant fraction of their less massive stellar components.They also show the lowest r h /R J and r t /R J ratios, implying that they are well inside their Jacoby radius.This would mean that their evolution is mainly dominated by their internal dynamics and not by the influence of the Galactic tidal potential.This result is in agreement with the N-body simulation's prediction that lobe Roche underfilling clusters will survive a number of relaxation times (e.g.Gieles & Baumgardt 2008).
NGC 6819 is a 1.9 Ga old and very massive system, M cl =19×10 3 M ⊙ , located at 8 pc from the Galactic centre with has lost less than 50% of its initial mass.NGC 2682 is older, 3.6 Ga, but less massive located at a similar height above the Galactic plane but about a kpc further, which have lost a somewhat larger fraction of its birth mass, ∼60%.NGC 2682 is slightly more concentrated, which may means that this cluster is more dynamically evolved, as evidenced by its dynamical age.In spite of their significantly different initial and nowadays masses, both systems have similar dissolution times due to the fact that the most massive one is also located at an innermost galactocentric distance.
On the other hand, NGC 2420, which almost the same age as NGC 6819, is located much further, R gc =10.7 kpc, and also at a larger distance from the plane, Z=869 pc, than NGC 2682 and NGC 6819.This cluster is twice and four times less massive than NGC 2682 and NGC 6819, respectively.However, its dissolution time is only a little lower mainly due to its location.These three clusters have similar r t /R J ratios of about 0.8, which means that they occupy almost all their Roche volumes, and therefore, their dynamical evolutions are modulated by the Galactic tidal field.
Finally, NGC 2425 is the furthest cluster in our sample but also the closest to the Galactic plane.It has an age similar to those of NGC 2420 and NGC 6819.This cluster shows clear hints of dissolution.On one hand, it has a r t /R J ratio larger than one, meaning that it completely fills its Roche volume since it may have been more exposed to external tidal forces (e.g.Ernst et al. 2015).It is also the system with the lowest concentration, which may mean that the stars that form this system are less gravitationally bound.In fact, we found a significant larger fraction of potential escapers for NGC 2425 than for the other clusters in our sample.Though, NGC 2425 is coeval, was formed with a similar initial masses, and it is located at similar galactocentric distance than NGC 2420, it has lost a significant larger fraction of its mass motivated by their different distances to the galactic plane.Taken into account its age, NGC 2425 will dissolve in only about 2.5 Ga from now.
Although our work is based in a small sample of only six clusters, we can obtain valuable conclusions which should be confirmed from larger samples.Our results reinforce the idea that the initial mass and the location in the Galaxy, mainly above or below the disc, play a fundamental role in the survival of open clusters.Those systems that do not fill completely their Roche volumes may survive longer, either they were born more concentrated or they lost a significant fraction of their outskirts.In fact, our results suggest that the concentration of the cluster increase with age.

Disc crossing times
In Sect. 3 we reported the existence of tails in several of the studied clusters, such as NGC 188.In order to investigate the influence of the Galactic plane in the morphology of the studied systems, we have estimated the last disc crossing.For this purpose, we use the orbits determined by Carrera et al. (2022) for all clusters except NGC 2425, for which we use Tarricq et al. (2021).In both cases, the orbits were derived using the PYthon galpy package (Bovy 2015) and the MW2014 potential.We refer the reader to the original papers for the details.We do not attempt to derive new orbits since the available values are the same, within the uncertainties, to those used in the orbit determination.
The times of the last disc crossing of each cluster are listed in the last column of Table 6.NGC 6819 has crossed the disc very recently, only about 12 Ma ago, which may be the explanation of its elongation.Moreover, it has relatively large disc pass frequency since it crosses the disc about every 50 Ma.Among the other clusters, NGC 6791 has passed the disc about 26 Ma ago, but there is no sign of elongation.The rest of the clusters have crossed the plane more than 30 Ma ago.In fact, the larger value is obtained for NGC 2425 but since this cluster do not reach a significant distance to the plane, its orbit is embedded inside the disc, and therefore, the last disc crossing is not relevant.
In order to evaluate this impact, we have computed the slopes of the radial density profiles of the studied clusters simply using a linear least square fitting.The results are shown in Fig. 6.Although the uncertainties in these external regions are large due to the small number of stars sampled, there is a hint of a correlation between the slope of the external region and the time that makes the clusters cross the disk for the last time without any clear relation with the actual distance to the Galactic plane.However, due to the small number of systems in our sample, this result should be confirmed by larger samples.

A cusp-core dichotomy for open clusters?
There are also clear differences in radial density profiles of the central regions among the studied clusters.In principle, these central regions are less affected by the external perturbations, and therefore, they may be dominated by the gravitational encounters of their members (e.g.O' Leary et al. 2014).A direct consequence of two-body relaxation is the redistribution of the stars, the mass segregation, and modification of their velocities due to the exchange of energy and angular momentum.In fact, the ages of the studied clusters are at least twice their t rh and therefore, the two-body relaxation has had enough time to act.All the clusters show clear signs of mass segregation, particularly in their cores, with mass segregation rations, Λ MSR , larger than two.The exception are NGC 188 and NGC 6791, with the lowest mass segregation ratios, which, as discussed above, should be related to the fact that these systems have lost a significant fraction of their less massive members.
Interestingly, the clusters with the lower Λ MSR value show flat profiles in their central regions, particularly clear in the case of NGC 6791.Meanwhile, for systems with higher values a power-law increase is apparent, the so-called cusp profile, which is specially clear in the case of NGC 2682.From our limited sample, the systems with a larger dynamical age show a core morphology.
The cusp-core dichotomy in the central density profiles is well-known among the globular clusters, more massive and older than the open ones (e.g.McLaughlin & van der Marel 2005).Although there are still some caveats about the mechanisms responsible for the cuspy profiles, it is widely assumed that clus-ters with such profile have had complex internal dynamics (e.g.Meylan & Heggie 1997).They are associated with systems that have suffered a gravothermal collapse, in which the increase of the central density becomes dramatic, the so-called core-collapse (Lynden-Bell & Wood 1968).However, the timescale for this event to happen is lower than in the case of globular clusters, so the oldest open systems could have experienced it (Spitzer 1987).
For globular clusters, several features have been related to cusp systems in comparison with core ones, such as a higher central density, a large binary fraction, a significant blue straggler population, and a high number of X-ray sources (e.g.Meylan & Heggie 1997).Binary fraction between 20% and 40% have been reported in the literature for the clusters in our sample (e.g.Bedin et al. 2008;Milliman et al. 2014;Thompson et al. 2021) but it becomes at 70%±17% for the central region of NGC 2682 (Geller et al. 2021) being 32%±3% for NGC 6791 (Bedin et al. 2008).However, Rain et al. (2021) reported the largest blue straggler populations for NGC 6791 and NGC 188 with 48 and 22 objects, respectively.NGC 2420, NGC 2682 and NGC 6819 contain 3, 11, and 15 blue straggler stars, respectively (Rain et al. 2021).These numbers are in agreement with our results, in spite our procedure is not the more adequate for these stars.From X-ray observations, it has been reported in the literature that NGC 2682 contains between 2 and 7 times more active binaries, normalised to the cluster mass, than NGC 6791, and even more with respect to NGC 6819 (van den Berg et al. 2013).
On the other hand, cored profiles can be a sign of the existence of massive remnants in the cluster centres, such as stellar mass black holes (e.g Merritt et al. 2004).Their presence is able to prevent the segregation of mass of the host systems in their central regions (e.g.Peuten et al. 2016)

Summary and conclusions
In this paper, we determine membership probabilities from Gaia astrometry for stars in the field of view of six of the oldest open clusters in the Milky Way in order to investigate their longevity: NGC 188, NGC 2420, NGC 2425, NGC 2682, NGC 6791, and NGC 6819.
From the study of their spatial distributions, we find that NGC 188 shows signs of a tail.NGC 2425 exhibits a nonisotropic distribution of the external region.In the case of NGC 6819, we notice an elongation of its central region.These structures are aligned with the directions of motion of each cluster.
The derived radial density profiles show some hint of flattening in the outskirts for two of the clusters, NGC 188 and NGC 2420.For the other, we only observe a change in the slope.There are also significant differences in the central regions.While NGC 2682 shows a power-law density increase, the so-called cusp profile, NGC 6791 exhibits a flat one, known as core profile.
We use LIMEPY and SPES set of models to characterise the observed radial density profiles.A fraction of potential escapers, stars with enough energy to leave the system or already out of it, is needed to properly reproduce the external regions.These models are used to determine some properties of the clusters such as their current masses, or half-mass, tidal and core radii.We found that the studied clusters are more extended than previously reported in the literature.
We discuss the obtained results, taking also into account the influence of the Galaxy.A high initial mass is needed to survive but also the location in the Galaxy, and mainly being at a large distance from the Galactic plane, play a role in the clusters' longevity.The most dynamical evolved systems do not fill completely their Roche volumes and lower r t /r c ratio.Moreover, all the cluster are segregated in mass but the most dynamically evolved clusters show a lower mass segregation ratio, which should be explained by the fact that these systems have lost a significant fraction of their less massive members.
Finally, the observed cusp-core dichotomy in the central regions may be related to different dynamical evolution.As in the case of globular clusters, cusp profiles may be related with more complex dynamical evolutions with a large presence of binaries and X-ray emission.On the other side, the presence of massive remnants, such as black holes, would prevent the mass segregation and formation of cusp profiles.More data, mainly kinematics, is needed to check these hypotheses.This will provide by Gaia for the brightest targets and by the massive spectroscopic surveys for the fainter ones.

Fig. 2 .
Fig. 2. Spatial distribution in Cartesian coordinates (x, y) of the studied cluster.The size of the points is related to the astrometric membership probability.The different populations are colour coded as in Fig. 1 excluding the binaries and BSS for clarity.The dark gray arrows in the bottom left corner show the direction of motion and are proportional to the velocity in each axis.The light gray arrows in the top right corner show the direction to the Galactic centre.
a) Fraction of potential escapers.tended than others.As a conclusion, King models are a good first approximation to the open clusters density profiles, especially if the external regions are not included.

Fig. 4 .
Fig. 4. Cumulative projected radial distribution of the different populations identified in each of the clusters, normalised to the total number of stars in each population.Error bars have not been plotted for clarity.

Fig. 5 .
Fig. 5. r h /R J (top), r t /R J (middle), and r t /r c (bottom) as a function of R gc , Z, τ rh , M cl /M ini , and cluster masses.Clusters studied here are labelled with different colours, as a function of ages, and symbols to facilitate their identification.Grey dots are the clusters studied by Angelo et al. (2023).

Fig. 6 .
Fig. 6.The slope of the outermost region of the radial density profile as a function of the last time the clusters crossed the disc.The clusters are colour coded as a function of their actual distance to the Galactic plane, using the same symbols as in previous figures.
. Initially, the formation of black holes in open clusters was discarded because of their low densities.However, recent dynamical studies have demonstrated that binary black holes can form in open clusters by different mechanisms (e.g Mapelli 2016; Kumamoto et al. 2020).In fact, very recently, it has been suggested by Torniamenti et al. (2023) the existence of a black-hole population to properly reproduce the density profile of the Hyades open cluster.

Table 1 .
Properties of the sampled clusters.

Table 3 .
Best-fitting parameters of LIMEPY models

Table 4 .
Best-fitting parameters of SPES models

Table 5 .
Summary of cumulative projected radial distribution of different populations.

Table 6 .
Physical properties obtained using clusters' masses derived from SPES models.