Stellar escapers from M67 can reach solar-like Galactic orbits

TL;DR

This study uses N-body simulations to explore if stars escaping from M67-like clusters can reach solar-like orbits, suggesting the Sun might have originated from such a cluster in the Milky Way.

Contribution

It provides the first detailed dynamical modeling of M67-like clusters' stellar escape patterns and their potential to produce solar-like orbit stars.

Findings

Escapers from all cluster types can have solar-like kinematics.

Depleted clusters produce the highest fraction of solar-like escapers.

Destroyed clusters have roughly twice as many solar-like escapers as surviving ones.

Abstract

We investigate the possibility that the Sun could have been born in M67 by carrying out $N$-body simulations of M67-like clusters in a time-varying Galactic environment, and following the galactic orbits of stars that escape from them. We find that model clusters that occupy similar orbits to M67 today can be divided up into three groups. Hot clusters are born with a high initial $z$-velocity, depleted clusters are born on cold orbits but are destroyed by GMC encounters in the Galactic disc, and scattered clusters are born on cold orbits and survive with more than 1000 stars at an age of 4.6 Gyr. We find that all cluster models in all three cluster groups have stellar escapers that are kinematicaly similar to the Sun. Hot clusters having the lowest such fraction $f_{\odot} = 0.06$ %, whilst depleted clusters have the highest fraction, $f_{\odot} = 6.61$ %. We calculate that clusters that are destroyed in the Galactic disc have a specific frequency of escapers that end up on solar-like orbits that is $\sim$ 2 times that of escapers from clusters that survive their journey.

Tables

Table 1. Table 1: Observational data and their transformation into Milky Way model parameters for the Sun and M67. References are: (a) Bellini et al. ( 2010 ) , (b) Schönrich et al. ( 2010 ) , (c) Reed ( 2006 ) and (d) Brunthaler et al. ( 2011 ) .

Observational parameter	Value	Reference
	M67
${\mu }_{\alpha }$ cos($\delta$)	-9.6$\pm$1.1 [mas yr-1]	a
${\mu }_{\delta }$	-3.7$\pm$0.8 [mas yr-1]	a
Radial velocity	33.78$\pm$0.18 [ km s-1]	a
Distance	815$\pm$81.5 [pc]	a
	Sun
$U_{\odot }$	-11.1$\pm$1.2 [ km s-1]	b
$V_{\odot }$	12.24$\pm$2.1 [ km s-1]	b
$W_{\odot }$	7.25$\pm$0.6 [ km s-1]	b
$z_{\odot }$	20$\pm$5 [pc]	c
$\delta R$	0.23 [kpc]	d

Table 2. Table 2: Results for our 16 M67 candidates. Coulmn list: number of stars remaining in the cluster, half-mass radius, mass, fraction of escapers that end up on solar-like orbits, fraction of escapers that end up on solar-like orbits that have escaped the cluster before the cluster age was 250 Myr, initial z 𝑧 z -velocity and guiding radius.

Cluster	$N_{\star }$	$r_h$	$m$	$f_{\odot }$	$f_{\odot ,250}$	$v_{z,0}$	$R_{g,0}$
		$\left[\mathrm{pc}\right]$	$\left[M_{\odot }\right]$	[%]	[%]	$[$ km s-1$]$	$\left[\mathrm{kpc}\right]$
A	2618	3.5	1707.0	1.67	0.38	-5.5	8.0
B	4208	3.3	2592.0	0.25	0.09	-12.3	7.4
C	4683	3.6	2812.4	1.07	0.61	4.4	8.1
D	183	2.8	174.7	6.61	1.03	-3.6	10.0
E	3209	3.3	2092.7	0.31	0.12	17.4	8.4
F	3361	3.4	2172.3	0.62	0.13	-8.8	7.9
G	6784	3.9	3757.4	0.49	0.22	14.7	7.4
H	9837	4.6	4994.8	0.38	0.13	-11.3	8.3
I	774	2.8	622.1	3.74	0.47	-5.8	8.0
J	5456	3.5	3184.5	0.06	0.02	20.7	8.2
K	639	2.17	545.2	2.02	0.24	17.4	7.2
L	1542	2.9	1141.4	0.77	0.17	-2.7	7.4
M	1277	2.6	980.2	5.07	0.27	1.5	9.0
N	379	2.4	980.2	2.84	0.47	2.1	7.7
O	256	2.3	244.4	2.10	0.64	-1.1	7.3
P	2951	3.5	1933.1	1.07	0.10	-18.7	8.4

Table 3. Table 3: The velocity dispersion of stars in the three galactic rectangular coordinates U , V 𝑈 𝑉 U,V and W 𝑊 W at a time of 4.6 Gyr. A comparison of our Milky Way model with the simulations of Gustafsson et al. ( 2016 ) and observations of Holmberg et al. ( 2009 ) can be seen. The ROSO simulation are similar to the Milky Way model, but where we use the GMC distribution from Rosolowsky ( 2005 ) which is more bottom heavy with a lower average GMC mass. ♠ ♠ \spadesuit : This value was only calculated for stars crossing the Galactic plane.

	${\sigma }_U$	${\sigma }_V$	${\sigma }_W$	$f_{400}$
	[ km s-1]	[ km s-1]	[ km s-1]	[%]
MW model	31.9	18.0	11.7	3.7
ROSO	29.9	15.7	10.5	2.1
Gustafsson et al. (2016)	37.5	18.9	16.9♠	1.8
Holmberg et al. (2009)	30.9	21.7	15.2	-

Equations

$$\rho (a)={\rho }_0(1-a),$$

$${\rho }_0(R)={\rho }_{02}{e}^{(-(R-R_s)/R_L)},$$

$${\rho }_{02}=\frac{3M_s}{2\pi a_0^2{c}_0{\sum }_{j=1}^N{e}^{-(R_j-R_s)/R_L}},$$

$$\frac{\mathrm{d}N}{\mathrm{d}{M}_i}\propto M_i^{-1.8}\mathrm{with}5.0\le {\log }_{10}(M_i/M_{\odot })\le 7.0,$$

$$M_t=\left[-0.25\cdot {\left(\frac{t-t_0}{10^7}\right)}^2+\frac{t-t_0}{10^7}\right]\cdot M_i,$$

$$R_c=20\cdot {\left(\frac{M_t}{5\times 10^5{M}_{\odot }}\right)}^{1/2}\mathrm{pc}.$$

$$d=100\cdot {\left(\frac{M_i}{10^6{M}_{\odot }}\right)}^{1/2}\mathrm{pc}.$$

$${\delta }_E=\frac{8G{M}^2{r}_h^3}{3m{b}^4{V}^2},$$

$$T_{ij}=-\frac{{\partial }^2\Phi }{\partial x_i\partial x_j},$$

Full text

††pubyear: 2019††pagerange: Stellar escapers from M67 can reach solar-like Galactic orbits–Stellar escapers from M67 can reach solar-like Galactic orbits

1 Introduction

M67 is an old, metal-rich open stellar cluster located at a distance of around 800 - 900 pc from the Sun. It is at a height of $z$ $\sim$ 400 - 450 pc (Sarajedini et al. (2009), VandenBerg & Stetson (2004), Friel (1995)) above the Galactic Plane which has allowed the cluster to survive for much longer than open clusters usually do, with an estimated age of 3.5 - 4.8 Gyr (Yadav et al., 2008).

The links between the Sun and M67 are numerous. M67’s estimated age is comparable to the age of the Sun, together with the cluster having a solar-similar chemical composition with [Fe/H] = 0.00 $\pm$ 0.06 (Heiter et al., 2014). A dominant fraction of stars analysed in M67 by Önehag et al. (2014) shows abundance patterns that are in better agreement with the Sun than solar twins in the Solar neighbourhood.

There are several indications from our Solar system that suggest a Solar birth cluster of substantial size such as M67 is needed. Meteoritic evidence of short-lived radioactive isotopes implies that the early Solar System has been enriched by an external source. A popular external source for enrichment is a supernova explosion with an preferred initial starting mass of $\sim$ 25 $M_{\odot}$ which has to be within a distance of $\sim$ 0.1 - 0.3 pc in order to provide enough enrichment for the Solar nebula (Adams, 2010). From cluster enrichment simulations with a total number of 2100 stars Parker et al. (2014) found that this enrichment process is able to enrich $\sim 1$ % of G-dwarfs that have always been single stars and that this process is quite stochastic.

The outer objects of the Solar System also provide clues to the Sun’s birth enviroment. With an large eccentricity of $\sim$ 0.84, Sedna has an unusual orbit compared to other solar system objects. Numerical studies (Morbidelli & Levison, 2004) have shown that Sedna’s orbital elements can be explained by a close encounter with a passing star where an impact parameter of $\sim$ 400 AU is needed.

The connection between M67’s position high above the Galactic plane and the Sun’s position has been addressed by Pichardo et al. (2012) who investigated their dynamical link by integrating the current positions of M67 backwards in time in a galactic model which included spiral arms and a Galactic bar. They argued that the Sun could not be born in M67, since the outer part of the solar system would have been destroyed as the Sun was ejected from the cluster onto the current solar orbit. This conclusion is based on the assumption that M67 was born on its current orbit. However, as shown by Gustafsson et al. (2016), if M67 was born in an orbit close to the Galactic plane and was then subsequently scattered to high altitudes by encounters with giant molecular clouds (GMCs), the kick needed for the Sun to escape M67 would not need to be so strong as to destroy the outer Solar System.

In this paper we investigate the orbits of stellar cluster escapers and the possibility that the Sun has been born in M67 by using a Milky Way model and $N$-body simulations of M67. We test the hypothesis that the Sun was tidally stripped away from M67 by interactions with the spiral arms and GMCs of the Milky Way as the cluster departed from its birth place in the Galactic plane to a height of $z$ $\sim$ 400 - 450 pc (Sarajedini et al. (2009), VandenBerg & Stetson (2004), Friel (1995)). The tidal stripping is caused by the galactic potential, which in our model contains both a static component and time-varying components that represent the spiral arms, Galactic bar and GMCs. Interactions with this time-varying force field can scatter clusters such as M67 up to high $|z|$ as shown by Gustafsson et al. (2016).

The simulations we preform can be divided into three steps. Step 1: We follow test particles, representing clusters, in the Milky Way model for 4.6 Gyrs and identify those that have similiar orbits to the current orbit of M67. We call these M67 candidates. Step 2: The evolution of each M67 candidate is computed in a $N$-body simulation where the effects of the galactic tidal field and encounters of the cluster with the GMCs are taken into account. Step 3: We follow the escaped stars from the $N$-body simulations in the Milky Way model til the present day and analyse their orbits.

The paper is organized as follows: In Section 2 we present the galactic model we use to represent the Milky Way and follow galactic orbits. In Section 3 we define M67-like candidates in terms of their orbital parameters and survivability based on their GMC encounters. Section 4 describes the parameters of the $N$-body model that is used to represent an initial model of M67. We present the tidal tensor which represents the galactic tidal field and explain how GMC encounters are incorporated. We show our results in Section 5 and discuss the different M67 model histories in terms of their orbital journey to the present day position of M67 and what orbits their stellar escapers end up on. We further discuss the results in Section 6 and give our conclusions in Section 7.

2 Milky Way Model

The Milky Way model follows Gustafsson et al. (2016) and consists of a axisymmetric potential (Binney, 2012), two spiral arms (Pichardo et al., 2003), a Galactic bar (Pichardo et al., 2012) and a distribution of GMCs along the spiral arms (Gustafsson et al., 2016).

The axisymmetric galactic potential in our model is based on Potential I of Binney (2012) which consists of a thin and thick stellar disc, a gas disc, and a stellar and dark matter spheroid representing the Bulge and the Halo. From the potential we removed the stellar spheroid representing the Bulge and replaced it with our model of the Galactic bar (see below). This gives a circular speed of $\sim$ 220  km s*-1*  at $R_{0}=8$ kpc.

2.1 Spiral arms

The spiral pattern consists of two spiral arms with a pitch angle of 15.5*∘* described by Pichardo et al. (2003). Each arm is represented by 100 inhomogenous oblate spheroids with semi-major and minor axes of 1000 pc and 500 pc, respectively. The spiral pattern starts at the radius where the Galactic bar ends at $R_{s}$ = 3.13 kpc with a distance between each spheroid of 500 pc. The spheroid density follows a linear density law

[TABLE]

where $a=(x^{2}+y^{2}+z^{2}/c_{0}^{2})^{1/2}$ and $c_{0}=500$ pc. This density law gives each spheroid a zero density at its boundary where $\rho_{0}$ is the central density. The central density for each spheroid has a exponential decline of

[TABLE]

with a radial scale height, $R_{L}$, of 3.9 kpc. $\rho_{02}$ is given in Pichardo et al. (2003) by

[TABLE]

where $M_{s}$ is the total mass in the spiral arms, $R_{j}$ is the Galactocentric distance of the spheroid’s centers and the sum is over one spiral arm. In our model we use a total spiral mass of $4\times 10^{9}$ $M_{\odot}$ and a constant pattern speed of 24  km s*-1*  kpc*-1*.

The force from each spheroid was pre-calculated in a grid as a function of radius $r$ and angle $\beta$ in the $x-z$ plane of the spheriod and used in a bicubic interpolation for the force calculation for each oblate spheroid representing the spiral arms. Due to symmetry $\beta$ was defined between $0^{\circ}$ and $90^{\circ}$ with a step size of $\Delta\beta=0.5^{\circ}$. The radius $r$ was defined between 0 and 20 kpc with a step size of 0.1 kpc, except for the inner 2 kpc where the step size was 0.01 kpc in order to get higher resolution for a more precise bicubic interpolation.

2.2 Galactic bar

The Galactic bar is represented by a prolate inhomogenous spheroid (Pichardo et al., 2003, 2012) with a linear density law similar to Eq. 1 with zero density at the boundary. The semi-major and minor axes of the Galactic bar is $a_{0}=3.13$ kpc and $c_{0}=1.0$ kpc which gives $a=(x^{2}/a_{0}^{2}+y^{2}+z^{2})^{1/2}$. The force exerted by the Galactic bar was pre-calculated in a grid as a function of $r$ and angle $\beta$. The radius in the grid was defined between 0 to 12 kpc with a step size of 0.1 kpc. The range of $\beta$ was between $0^{\circ}$ and $90^{\circ}$ with a step size of $\Delta\beta=1.0^{\circ}$. The numerical integration which was used to produce the pre-calculated grid was done in a similar fashion as for the spiral arms. The total mass of the Galactic bar was set to $1.6\times 10^{10}$ $M_{\odot}$ with a constant pattern speed of 55  km s*-1*  kpc*-1*.

2.3 Giant Molecular Clouds

The implementation of GMCs follows Gustafsson et al. (2016) where the maximum mass of a GMC, $M_{i}$, was selected according to the distribution function

[TABLE]

based on the observations and simulations of Williams & McKee (1997) and Hopkins et al. (2012), respectively. We give each GMC a lifetime of 40 Myr corresponding to a few free fall times. During its lifetime the GMC will increase its mass to a value of $M_{i}$ over 20 Myr and then decrease back to zero over the next 20 Myr with a parabolic mass evolution of

[TABLE]

where $t_{0}$ is the birth time of the GMC and $t$ is the current time, both in years. The initial birth times for the first GMCs were chosen randomly between the time -40 and 0 Myr. When the life of one GMC ended a new is born such that the number of GMCs in our simulations is always constant. We adopt the total Galactic molecular mass of $1.0\times 10^{9}$ $M_{\odot}$ from Williams & McKee (1997) resulting in a constant number of 2500 GMCs. Each GMC is represented by a Plummer (1911) sphere with a core radius of

[TABLE]

Newborn GMCs are distributed within a distance of $\pm$ 50 pc from the median line of the spiral arms in the Galactic plane, at Galactocentric distances 4 kpc $\leq R_{0}\leq$ 9kpc, with a number density decreasing exponentially along the arms similar to Eq. 2 in order to match the density decrease of the arms. Each GMC initially moves with the local circular rotation velocity plus a velocity scatter which has a Gaussian distribution of $\sigma=7$  km s*-1*  for each of the three velocity components. The total surface density profile for the Milky Way model can be seen in Figure 1. The position of the Sun is represented by the yellow star and the black dots are the GMCs at a time of 4600 Myr into the simulation, corresponding with our present time.

2.4 Matching Milky Way model with observations

The test particles representing our clusters were randomly placed in the Galactic disc with Galactocentric distances, 4 kpc < $R$ < 10 kpc. Like the GMCs, the particles were given a gaussian velocity distribution of 7  km s*-1*  in each of the three velocity directions, together with an initial local circular rotation speed. The particles were followed for 4.6 Gyr with a total number of test particles of 25000 divided into five independent runs. Each run was done in a fixed reference frame with a rotating Galactic bar and spiral arms. A constant timestep of 0.1 Myr was used for all Milky Way model runs.

Our model was able to reproduce the observed velocity dispersion of solar-type stars within 300 pc of the Sun by Holmberg et al. (2009) as a function of time seen in Figure 2. The test particles used for calculating the velocity dispersion were taken from within 300 pc of a Galactocentric radius of 8 kpc for all azimuthal angles in order to get better number statistics. We see that the evolution of $\sigma_{U}$ and $\sigma_{V}$ follows the observations quite well. The evolution of $\sigma_{W}$ does not follow the observations and is displaced by a few  km s*-1*  which suggest that our Milky Way model is biased against test particles going onto high altitudes. For the purpose of our investigation, this means that scattering test particles into M67-orbits in our Milky Way model will be more difficult and thus will most likely underestimate the number of M67 candidates in our model.

2.5 Solar-like orbit

In order to compare the orbits of escaped stars with that of the Sun we first need to define a Solar-like orbit. We use three parameters to classify orbits in the Milky Way model: the eccentricity $e$, the maximum height above the Galactic plane $|z|_{\mathrm{max}}$ and the guiding radius $R_{g}$. The eccentricity was calculated as $e=(r_{a}-r_{p})/(r_{a}+r_{p})$, where $r_{a}$ and $r_{p}$ are the apoapsis and periapsis of a single orbit, which we define as a full rotation of 360*∘* in the Milky Way disc. $|z|_{\rm{max}}$ is defined as the maximum $|z|$ of a test particle during a full orbit. The guiding radius was defined as $R_{g}=L_{z}/v_{c}$, where $L_{z}$ is the $z$-component of the angular momentum and $v_{c}$ is the local circular velocity at the test particle’s current position. To construct the orbit we used the local standard of rest and the current $z$-position of the Sun. We adopted the uncertainty of the Galactocentric distance $\delta R$ from Brunthaler et al. (2011) with a Galactocentric distance of 8 kpc as defined in our Milky Way model, the parameters of which can be seen in Table 1. We integrated the orbit of the Sun 300 Myr back in time to obtain $e$, $|z|_{\mathrm{max}}$ and $R_{g}$. No GMCs were used for the backwards integration since this could result in an unrealistic higher spread of the three parameters. The backwards integration of the Sun was done for 5000 test particles with initial values sampled according to a Gaussian distribution of the values in Table 1. From this we obtained $e_{\odot}=0.1$ and $|z|_{\mathrm{max},\odot}=98.7$. In order for cluster escapers to end up on solar-like orbits, they need to be dynamically colder than the current values of $e_{\odot}$ and $|z|_{\mathrm{max},\odot}$ and not necessarily match these values. We do not require the guiding radius to match the Solar value as a criteria for escapers to be classified as solar-like. This is because $R_{g}$ is mainly modified by churning of orbits by the spiral arms and not scattering off GMCs. The spiral arms are transient in the Milky Way but static in our model, which means that the position in terms of guiding radius of our escapers could be quite different depending on the spiral arm evolution of the Milky Way during the last 4.6 Gyrs. In regards to how normal it is to be solar-like in our model we compare the final position of our 25000 test particles to our definition of what we classify as a solar-like orbit. To the left in Figure 3 the position of all 25000 particles at a time of 4.6 Gyr are shown, with the red rectangle defining the area where test particles are solar-like. To the right in Figure 3 we only see the 8578 particles that have a guiding radius between 7 and 9 kpc. The fractions of particles that have solar-like orbits at the end of the simulations is quite high, 13 % in the full sample and 11 % in the subset with 7 $\leq R_{g}\leq$ 9 kpc. We would therefore expect that a high fraction of escapers, from clusters that are destroyed in their early life in the Galactic disc, will be on solar-like orbits at a time of 4.6 Gyr.

3 Cluster survival and M67-like orbits

We define M67-like orbits in the same way as we define the orbital parameters for the Sun. We constructed the orbit using the radial velocity, proper motion and position of M67 seen in table 1. With the absence of GMCs we integrated the orbit of M67 300 Myr back in time to get $e$, $|z|_{\mathrm{max}}$ and $R_{g}$. We found the orbital parameters of M67 to be $R_{g}=7.89\pm 0.27$ kpc, $|z_{max}|=480\pm 56$ pc and $e=0.12\pm 0.02$. Out of the 25000 test particles in our Milky Way model, 37 of them had similar orbital parameters as M67 to within one standard deviation and can be seen as green points in Figure 4 where the location of the initial cluster test particles can be seen in blue. The red points are the position of the cluster particles after 4600 Myr.

GMCs can help scatter clusters to high galactic orbits, but they can also destroy the clusters via tidal forces if the GMC comes too close. We therefore had to evaluate which of our 37 cluster particles would survive their journey to the location of M67 over a time span of 4.6 Gyr. We consider encounters between a GMC and a cluster particle to be close when the closest separation $r_{\mathrm{sep}}\leq d$ where

[TABLE]

Following Eq. (19) from Gustafsson et al. (2016), the fractional change in kinetic energy in a cluster as it encounters a GMC can be approximated as

[TABLE]

where $G$ is the gravitational constant, $M$ is the mass of the GMC, $r_{h}$ is the half-mass radius of the cluster, $m$ is the mass of the cluster, $V$ is the relative speed between the GMC and cluster and $b$ is the impact parameter. If the cumulative $\delta_{E}>1$ for all the GMC encounters, the cluster is tidally destroyed. For the half-mass radius of M67 we chose, in accordance with previous investigations of M67 in Hurley et al. (2005) a value of 4 pc and 36000 stars corresponding to a mass of $\sim 2.2\times 10^{4}$ $M_{\odot}$. The values of the half-mass radius and mass will change as the cluster evolves but remained constant in Eq. 8 as a conservative estimate of the cluster survivability. Out of the 37 M67-like clusters we estimate that 16 will survive their GMC encounters. We classify these 16 clusters as M67 candidates and further investigate their histories.

4 $N$-body simulations

In order to obtain realistic escape histories of our 16 M67 candidates we carried out $N$-body simulations for each cluster. The GMC encounter histories and orbits in the galactic field of our Milky Way model will shape the evolution of each M67 candidate differently. To follow the galactic tidal enviroment of each cluster we used the $N$-body code nbody6tt (Renaud et al., 2011). nbody6tt is a modified version of nbody6 (Aarseth, 2003) where a tidal tensor can be implemented to described the tidal forces the cluster experiences during its orbit.

4.1 Initial parameters of M67

For the initial parameters of M67 we adopt initial values similar to those used by Hurley et al. (2005), with 36,000 stars and a half-mass radius of 4 pc. The stars were distributed in mass according to a Kroupa (2001) IMF in a range of 0.1 - 50 $M_{\odot}$ and spatially distributed according to a Plummer (1911) distribution. Stellar evolution was turned on with a solar metallicity ($Z=0.02$).

4.2 Tidal tensor

The tidal tensor is defined in Renaud et al. (2011) as

[TABLE]

where $\Phi$ is the gravitational potential of the galactic enviroment at the cluster’s position. The tidal tensor was constructed by following each M67 candidate through the Milky Way model. We calculated $T_{ij}$ by numerical differentiation of the force for each time step within a cube of size $2\delta$ centered on the cluster particle. The value of $\delta$ was chosen to be 15 pc, corresponding to a few half-mass radii of the cluster, such that the tidal tensor would be able to resolve the fine structures of the gravitational potential of the Milky Way model, but still not too large such that it would introduce errors into the tidal tensor from GMCs coming too close.

4.3 Implementing GMCs into nbody6tt

The 16 M67 candidates that are expected to survive their GMC encounters have on average 43.6 registered GMC encounters according to Eq. 7. Based on the work done by Gustafsson et al. (2016) we know that GMC encounters that inject less than a few $\delta_{E}=0.01$ will not affect the cluster in a significant way. We therefore only focused on the GMC encounters with $\delta_{E}\geq 0.01$ when it came to implementing them into the $N$-body simulation. The GMC encounters with $\delta_{E}<0.01$ were still considered but because of their weaker interactions they are treated as part of the galactic tidal field defined by the tidal tensor. By considering only the stronger GMC encounters we see that our 16 clusters on average have 4.9 GMC encounters.

The GMCs in nbody6tt behave in a similar fashion as in the Milky Way model in terms of being a Plummer (1911) potential with a mass evolution and core radius according to Eq. 5 and 6, respectively. From the Milky Way model we record the minimum separation between the cluster and GMC during each encounter, $r_{\mathrm{sep}}$, the relative position $\vec{r}(r_{\mathrm{sep}})$, the relative velocity $\vec{v}(r_{\mathrm{sep}})$ and time of separation $t(r_{\mathrm{sep}})$. Using this information we integrate the relative orbit of the GMC and cluster back to the birth time of the GMC to obtain the relative initial position and velocity of the GMC. This integration was needed because the mass of the GMC changes as a function of time. A picture of the scenario of the backwards integration can be seen in Figure 5.

4.4 Escapers

The treatment of stellar escapers in nbody6 is normally done by demanding that the stars have a postive orbital energy and are located more than 2 tidal radii away from the cluster. Since we can have compressive tides in nbody6tt, the tidal radius is no longer a good criterion for our escapers (Renaud et al., 2011). The criteria we instead adopted for a star being an escaper was that it should have a postive orbital energy and be at least ten half-mass radii from the cluster.

As the $N$-body simulations of the M67 candidates evolve we record the escape time, $t_{\mathrm{esc}}$, and velocity, $\vec{v}_{\mathrm{esc}}$, for each stellar escaper. Each cluster was then re-run in the Milky Way model, where the orbit was integrated with the same initial values and GMCs as earlier. Hence the cluster follows the same orbital path as before. The stars escape the cluster as test particles following their escape history from the $N$-body simulation. As the stars escape they will still feel the gravitational force from the cluster, which is itself evolving in mass and half-mass radius. We represent the force from the cluster as a Plummer (1911) potential with a scale paramter equal to $a=r_{h}/1.305$, with a linear interpolation of the mass and half-mass radius as a function of time from the results of the $N$-body simulation.

5 Results

In this section we first focus on a specific model cluster, which we refer to as A, and then later give a more generalized view of all the clusters and trends we see. Cluster A was chosen because it resembles M67 today, and has one of the highest fractions of escapers that end up on solar-like orbits, $f_{\odot}$.

5.1 Cluster A

Cluster A’s orbital evolution consists of five GMC encounters, four of which occur before the cluster is 1 Gyr old. The first two encounters are only seperated in time by 2 Myr with encounter times of 518.8 and 520.8 Myr. Similarly the two next encounters also happen in a pair with only 0.9 Myr between them, at times of 819.8 and 820.7 Myr. These four GMC encounters cause the cluster to lose a significant number of stars which can be seen in Figure 6. Here the number of stars in the cluster can be seen, together with vertical shaded regions that represent the lifetime of each GMC for each GMC encounter. Roughly 40 % of the stars are in fact lost within the first 1 Gyr which means that the stellar escapers are leaving the cluster while it is still in a cold disc-like orbit. Between GMC encounters, e.g. in the range of $t=[1000,3300]$ Myr, we see a constant loss of stars that is a consequence of the internal dynamics of the cluster and the external galactic tidal field. In this time between the GMC encounters, the cluster loses around 5 stars per Myr.

As the cluster experiences GMC encounters and the galactic tidal field, we see in Figure 7 that after the initial expansion, typical of $N$-body models with mass loss, the half-mass radius is fairly constant at a value around 5 pc, from a couple of hundred Myr and untill the cluster is around 3 Gyrs old. The spikes we see in this time span are produced by the GMC encounters which causes the half-mass radius to increase as a consequence of stars being tidally stripped from the cluster. Ultimately these stars escape and the half-mass radius settles down to a value of $\sim 5.5$ pc. After the last GMC encounter in Figure 6 and 7 we still see a constant loss of stars, however the half-mass radius continues to drops down to $\sim 3.8$ pc before stabalizing.

Because stars escape from cluster A at very different times their orbits evolve separately after they are lost from the cluster. Figure 8 shows the orbits of stars that escape as a consequence of three GMC encounters for cluster A. The left panels show the stars at their escape time and right panels show the orbits of the same stars at the end of the simulation. From the left panels we see how the initial $|z|_{\mathrm{max}}$ of the escapers depends on when the stars escape. This is a consequence of the orbital heating the cluster undergoes as its orbit evolves in the Galactic enviroment and is scattered to higher $|z|$. The initial orbit of an escaped star is therefore very dependent on what kind of orbit the cluster is on. This is because a majority of the stellar escape speeds relative to the cluster are below 5  km s*-1*  which we will discuss in later sections.

From the two first rows of panels, for stars that escape before 1 Gyr, we see that the initial structure of the escapers in $e$-$|z|_{\rm{max}}$ space is completely washed out by the present day due to their interactions with the galactic enviroment. We therefore see no distinction between the final orbits because the stars have been lost early enough in the cluster’s history to be equally affected by the galactic enviroment. For the last encounter that happens at $t=3423$ Myr we see that the escapers have not had enough time to fully diffuse out from the cluster’s original orbit. These stars also escape the cluster at a time where the cluster is already at a high $|z|$ orbit which means that the rate of interactions with spiral arms and GMCs are significantly decreased, causing the diffusion of escapers to be slower.

Figure 9 shows the present day orbital parameters of all the escapers from Cluster A. The white star represents the current orbit of M67 and the red rectangle indicates the area for which we defined the stars to be on a solar-like orbit. We see a large concentration of stars around the cluster, which indicates that the stars that have escaped recently, still have orbits similar to the cluster’s. This makes sense, since at later times the clusterwill be at high altitude and therefore the interaction with the Milky Way disc will be less frequent and the diffusion of the stars from the cluster will take longer. About 1.67 % of stars that escape from cluster A end up on solar-like orbits, the second largest fraction for the clusters that still have a reasonable number of stars left. The reason for cluster A’s numerous solar-like escapers can be found when investigating when the cluster experiences its GMC encounters.

Looking at Figure 10 we see the fraction of escapers that end up on solar-like orbits, $f_{\odot}$. Here $f_{\odot}$ is shown as function of time with a timestep of 100 Myr. We see that $f_{\odot}$ is highest in the begining of the simulation where the cluster is still in a cold disc-like orbit seen in Figure 11 which shows the $z$-position of the cluster as a function of time. The GMC encounters do not change this fraction, however they substantially cause a lot of stars to escape during the time while the cluster is still in this cold disc-like orbit. As the cluster is gradually heated by interactions with the spiral arms and GMCs, the fraction of solar-like orbits decreases. After around 2 Gyr the cluster is already at the same height as M67 and we see that the fraction of ejected stars that eventually reach solar-like orbits is non-existent after this point. An important factor for a M67-like cluster to have a high $f_{\odot}$ value seems to be that it has early GMC encounters while the cluster is still in a cold disc-like orbit.

Aside from having a roughly constant stellar escape rate between GMC encounters, Figure 12 shows signs of a periodic signal in the escape rate. Here we see the number of escapers as a function of their escape time which shows a periodic signal with a amplitude of $\sim 250$ stars and a period of $\sim 800-900$ Myr. This period agrees well with the timescale of closest approach between the cluster and the spiral arms at the current $R_{g}$ of the cluster. This signal could therefore be a manifestation of the two spiral arms in the galactic model, however this is hard to verify since the spiral arms are stretched over a large area. From Figure 12 we also see GMC encounters lead directly to the loss of around 10 % of the cluster’s stars per encounter. Even so, the majority of stars are not lost as a direct result of an GMC encounters but instead to the internal dynamics of the cluster and its interaction with the galactic tidal field. Approximately 80 % of the stellar population is lost this way.

5.2 Comparison of three cluster histories

In order to explore the spread of cluster histories in our sample we compare the evolution of three clusters: A, B and C. The three cluster models are chosen based on their high values of $f_{\odot}$, their different GMC encounter histories, and the fact that all three clusters survive. Even though they have quite different evolution in terms of GMC encounters the clusters still end up with roughly the same size and number of stars after 4.6 Gyr. The values of $f_{\odot}$ for the clusters are, however, considerably different which can be attributed to the different timings of these GMC encounters. Results for all 16 M67 candidates can be seen in Table 2.

Cluster B has one very strong GMC encounter which causes the largest expansion of the half-mass radius of all three clusters to a value of around 8.2 pc seen in Figure 7. Afterwards the half-mass radius steadily decreases until the end of the simulation. When this single GMC encounter occurs the cluster has already been heated up to $|z|_{\mathrm{max}}\sim$ 300 pc as seen in Figure 11. The encounter causes around 8000 stars to escape, but since the cluster is already at a high $|z|_{\mathrm{max}}$ the number of stars that end up on solar-like orbits are considerably lower for cluster B than the two other clusters, with $f_{\odot}=0.25$ %. It should be noted that cluster B starts on a hotter orbit than A and C, with $v_{z,0}=-12.3$  km s*-1*  which corresponds to $|z|_{\mathrm{max}}\sim 200$ pc (see Figure 11). Cluster B has around the same escape rate as A and C, seen in Figure 12, but because B is born on this hot orbit the fraction of escapers onto solar-like orbits at $t_{\rm{esc}}=0$ Myr are a factor of 4 and 2 smaller than A and C, respectively. This can be seen in Figure 10.

Cluster C has two GMC encounters within the first 100 Myr of its evolution. This causes the cluster to lose around 3000 stars. At around 2400 Myr the cluster has a third encounter, however, this encounter seems to be particularly weak and we do not see an associated loss of stars escaping. A slight increase in the half-mass radius suggests a small expansion of the cluster due to the encounter. The next encounter happens at around 2900 Myr and like the former encounter it does not seem to have a great effect on the cluster. The last encounter happens just before the simulation ends and thus these escapers end up on similar orbits to the cluster. Just like clusters A and B, cluster C has roughly the same stellar escape rate while it starts out on a cold orbit like A with $v_{z,0}=4.4$  km s*-1*. The difference between cluster A and C, is that C does not have any early encounters which means that the $f_{\odot}$ in the early stages of the cluster history is lower than A, seen in Figure 10. Since this is the most important part of the cluster’s evolution in terms of losing stars that end up on solar-like orbits, $f_{\odot}$ is only 1.07 % for cluster C.

When looking at Figure 12 we see quite different evolutions. As mentioned before, we see a possible indication of the spiral arms affecting cluster A, however, it is not that clear for cluster B and C. The reason for this might be a result of the difference in the evolution of the eccentricity of the cluster’s orbits. Cluster A evolves with a fairly stable eccentricity of 0.08 and thus the spiral arms have a more periodic effect on the cluster. Cluster B has a eccentricity of about 0.14 for its first 1800 Myr until it has its GMC encounter which lowers the eccentricty to $\sim 0.11$ for most of the remaining simulation. Cluster B starts out not showing any periodic signal, but soon after the GMC encounter it starts to show a signal somewhat similar to cluster A even though it appears to be significantly weaker. Cluster C experiences an eccentricity of $\sim 0.15$ for most of its orbit and thus the signals from the spiral arms are scrambled together. The average stellar escaper rates of cluster A, B and C are quite similar when looking at the escape rates between their GMC encounters. The majority of stars lost for all three clusters are a result of the galactic tidal field and the internal evolution which means they end up with relatively similar cluster parameters.

Figure 13 shows $R_{g}$ as a function of time for clusters A, B and C. Cluster A experiences the most churning with changes ranging from 8.7 to 7.4 kpc. Cluster A is therefore more affected by the spiral arms of the Milky Way model which was also our suspicion when looking at the periodic signal in Figure 12. For clusters A and B we see a small change in $R_{g}$ at the time of their GMC encounters where as it is not as clear from cluster C. We stated earlier that we did not consider $R_{g}$ when classifying what a solar-like orbit was in our Milky Way model. Figure 13 shows exactly why this would not make sense to do, since the churning of the cluster and stellar escapers depend on what spiral structure is present and if it is transient or not.

5.3 Trends across all 16 models

We now expand our view and look at all of the 16 M67 candidates. In Figure 14 the number of stars in the cluster after 4600 Myr is shown as a function of the initial $z$-velocity of the cluster. Each cluster is marked and colour coded in regards to the number of GMC encounters it experiences. The vertical line indicates the initial 7  km s*-1*  dispersion used in the Milky Way model. The vertical line marks the $N_{\star}=1000$ threshold below which we define M67 candidates to be depleted of stars. The threshold is roughly based on the observed luminous mass of M67 of 1270 $M_{\odot}$ (Fan et al., 1996). From Figure 14 we can divide the 16 M67 candidates up into three groups labeled with the following symbols:

• $\triangle$

Hot clusters - Clusters that are born with a high initial $z$-velocity.

• $\Diamond$

Depleted clusters - Clusters that are born on cold orbits but are destroyed by GMC encounters in the disc.

• $\circ$

Scattered clusters - Clusters that are born on cold orbits with $N_{\star}\geq 1000$ at an age of 4.6 Gyr.

Given that the initial $v_{z}$ distribution of our cluster particles was Gaussian with $\sigma_{z}=7$  km s*-1*, we can see that there is a bias towards clusters born with high $|v_{z,0}|$. Hot clusters get to high $z$ without having to go through as many GMC encounters and risk getting destroyed. Depleted clusters and scattered clusters are in a sense more interesting since they are not outliers in the initial $|v_{z,0}|$ distribution. These clusters all start out with a low $|v_{z,0}|$ and therefore need more heating in terms of GMC encounters to get to the present day $z$-position of M67. We see that the depleted clusters go through the most GMC encounters, between 8 and 10, and as a result they are depleted of stars to the point where they do not match M67. The scattered clusters are born on dynamically cold orbits and manage to survive their 3 - 5 GMC encounters.

Figure 15 shows the 16 M67 candidates with $f_{\odot}$ as a function of the time of their first GMC encounter, $t_{\mathrm{enc,first}}$. The colour coding and symbols are identical to those in Figure 14. We see that hot clusters all take longer than 900 Myr before they experience their first GMC encounter and have $f_{\odot}$ below $\sim 1$ %, except for cluster K which is an outlier with $f_{\odot}=2.02$ %. Cluster K experiences an early GMC encounter compared to the other hot clusters and also loses significantly more stars, which could explain the higher value of $f_{\odot}$. Depleted clusters and scattered clusters all experience their first GMC encounter before 1 Gyr and range in $f_{\odot}$ from 0.77 to 6.61 %. From Figure 14 and 15 there seems to be an anti-correlation between the initial $|v_{z,0}|$ and $f_{\odot}$. We also see that clusters that start with a high $|v_{z,0}|$ generally take longer before they experience their first GMC encounter because they spend most of their time away from the Galactic disc where GMC encounters occur.

The final position in the $e$-$|z|_{\rm{max}}$ plane for the escapers associated with clusters A to P can be seen in Figure 16. Depleted clusters are shown first, followed by scattered clusters and hot clusters. For the depleted clusters we see that the escaped stars are located closer to the Galactic disc where the stars are lost due to GMC encounters. For clusters D, O and N we see a lower concentration of escapers around the present day position of M67 indicating that they have lost a higher fraction of their stars before they settled on this orbit. The scattered clusters have their escapers less concentrated to the Galactic disc and a higher concentration around the present day position of M67. This is because they retain more of their stars on their journey towards the present day orbit of M67. The hot clusters have their concentration of escapers higher above the Galactic plane which result in them having a lower $f_{\odot}$ value. They have fewer GMC encounters and less interaction with the Galactic disc which means that the diffusion of stellar escapers is less pronounced for some of the cluster, namely H, G, E and J. Again cluster K, which experiences its first GMC encounter at $\sim 100$ Myr, is an outlier.

6 Discussion

In order to compare our Milky Way model to observations and other simulations we can look at Table 3. When comparing our results to the results of Holmberg et al. (2009) we see that $\sigma_{U}$ and $\sigma_{V}$ match the observations quite well, while $\sigma_{W}$ shows the largest discrepancy. As mentioned earlier this suggest that our Milky Way model is biased against test particles going onto high altitudes. For the purpose of our investigation, this means that scattering test particles into M67-orbits in our Milky Way model will be more difficult and thus will most likely underestimate the number of M67 candidates in our model. This does, however, also mean that stellar escapers will experience less vertical scatter which could increase the number of solar-like escapers we see. It should be noted that the results from Holmberg et al. (2009) which used the Hipparcos parallaxes to improve the accuracy of the Geneva-Copenhagen Survey (Nordström et al., 2004) is not complete to within a distance of 300 pc which we use to classify the Solar-neighbourhood. As stated in Nordström et al. (2004) the brightest F stars are complete to a maximum distance of $\sim 70$ pc and G5 stars to within $\sim 40$ pc.

In comparison, the fraction of test particles that are located at $|z|$ > 400 pc, $f_{400}$, at the end of the simulation is 3.6 % compared to 1.8 % from Gustafsson et al. (2016) which can be seen in Table 3. Both of these values have been calculated after taking the initial Galactocentric radius into account and weighing the different test particles following the exponential gas disc of the Milky Way model of $\sim\mathrm{exp}(R_{i}/4.8\,\mathrm{kpc})$ as the test particles where uniformly distributed in both simulations. These two values do not agree that well, however, in comparison the simulation of Gustafsson et al. (2016) suffer from small number statistics where 500 and 1000 test particles were used, whereas our simulations are based on five independent runs with a total of 25000 particles. The initial birth radii of the test particles will also have an effect on the discrepancy between our model and the model of Gustafsson et al. (2016). Their initial birth radii were defined between 4 and 9 kpc whereas ours are between 4 and 10 kpc resulting in a higher fraction of stars being able to reach higher $|z|$ due to the lower strengh of the Milky Way potential at this Galactocentric radius.

We also simulated a different galactic model (which we will refer to as ROSO from now on) with a different GMC distribution from Rosolowsky (2005) with $\propto M_{i}^{-1.5}$ making it less top heavy in the mass interval of $4.0\leq\mathrm{log}(M_{i}/M_{\odot})\leq 6.5$ resulting in a constant number of 8770 GMCs. All of the velocity dispersions are lower than the values from our standard the Milky Way model and the results from Holmberg et al. (2009). In particular $\sigma_{V}$ is lower since the lower mass GMCs are not able to increase the velocity dispersion to match observations. The lower average GMC mass results in a lower $f_{400}$ of 2.1 % compared to our Milky Way model.

Pichardo et al. (2012) have claimed, based on simulations done in a galactic potential very similar to our Milky Way model (but without GMCs), that the Sun could not have been born M67. The arguement presented by Pichardo et al. (2012) was that in order for the Sun to be kicked out of M67 to a very different solar orbit would damage the outer parts of the Solar system. In their simulations they integrated M67 backwards in time from its current location. The result was that the Sun would have to leave M67 with a velocity larger than 20  km s*-1*  in order for the Sun to end up on the current solar orbit, since this is roughly the $v_{z}$ by which M67 currently crosses the Galactic plane. Kicks produced by three-body encounters of this magnitude would destroy the outer Solar System and the probability of an GMC being responsible for ejecting the Sun from M67 was calculated to be smaller than $10^{-4}-10^{-5}$. The vast majority of stellar escapers in our simulations all escape with speeds lower than 5  km s*-1*  which can be seen in Figure 17 for cluster A. Here we have divided the escapers up into two groups representing stars that escape as a consequence of an GMC encounter and stars that escape between GMC encounters due to internal and tidal effects. We see that most of the stars are below 10  km s*-1*  with the majority being around a few  km s*-1*  for both distributions. We see that GMC encounters cause a wider distribution, but drops off rapidly at $\sim$ 10  km s*-1*  unlike the stars that escape between GMC encounters. These low escape velocities decreases the risk of destruction of the Solar System compared to the results found by Pichardo et al. (2012) which hardly had any escape speeds below 20  km s*-1*  which would be enough to cause this destruction. This is because we assume that the cluster is born in the Galactic disc and via interactions with the galactic field and GMCs M67 is scattered to its current orbit and not initially born in it while our escaper orbits also diffuse out from the orbit of the cluster. The relative velocity of around 20  km s*-1*  is therefore not needed in our scenario.

To investigate the effect of GMC encounters on the Solar System we examine the tidal work from GMC encounters by using Eq. 8 to calculate the fractional change in energy of the orbit of Neptune during a GMC encounter. We pick the most violent GMC encounter occuring for the 16 M67 candidates with $M_{t}=4.3\times 10^{6}M_{\odot}$, $r_{\rm{sep}}=25.8$ pc and $v(r_{\rm{sep}})=43.7$  km s*-1*, and adopt Neptune’s apocenter of 30.33 AU to maximize the potential destruction of the orbit. We find $\delta_{E}=8\times 10^{-8}$ and can therefore conclude that GMC encounters can shape the internal and orbital evolution of stellar clusters, but they will have an negligible effect on the Solar System.

There is also a time limit when it comes to how long the Sun should be able to stay in a cluster like M67. With a stellar density of $\sim 1000$ pc*-3*, which is significantly higher than our representation of M67, Adams (2010) showed that the Sun would have to leave M67 before the cluster has reached an age of 250 Myr in order not to perturb the orbits of Uranus and Neptune by stellar flybys. Even if we might have a larger time limit, we can use the time limit of 250 Myr as a tight constraint and look at how $f_{\odot}$ changes if we only consider stars that escape the M67 candidates before a time of 250 Myr. This is represented by the parameter $f_{\odot,250}$ in Table 2. We see that clusters such as A, B, C, E and G all lose between 20 to 45 % of their stellar escapers that end up on solar-like orbits before the time limit of 250 Myr. The majority of our solar-like escapers, are lost early where the cluster is on a colder, disc-like orbit. Even restricting ourselves to stars escaping before 250 Myr we find escapers from our model clusters that reach solar-like orbits at the present day.

The number of GMC encounters for each cluster ranges from 10 down to zero encounters. Some of the clusters only have 0 - 2 encounters and even so, they still reach the height of M67. This is because they are born on dynamically hot orbits with initial $|v|_{z,0}>7$  km s*-1*. It is not clear whether or not it is in fact possible for a rich M67-like cluster to be born a high $|z|$ orbit, i.e. at the tail of the $z$-velocity distribution. We also see that clusters that are born on dynamically hot orbits have smaller $f_{\odot}$ and therefore make them less likely to be the birth place of the Sun.

We focus therefore on the two other groups of M67 candidates, the scattered clusters and depleted clusters. These are the ones which experience GMCs encounters and survive, i.e. they have a healthy number of stars at the end of the simulation, and the ones which have a low number of stars, effectively meaning that these clusters are destroyed in the disc before they reach a M67-like orbit. It should be noted that the distinction between depleted clusters and scattered clusters is approximate, since there is a continous variation of properties between the clusters in these two groups. From looking at Figure 14 we classify cluster A, C, L and M as scattered clusters and D, I, N and O as depleted clusters. Using this classification we calculate that it is $\sim$ 2 times more likely that a solar-like escaper will be produced by a cluster that gets destroyed in the Galactic disc before it gets scattered up to an M67-like orbit versus a cluster that survives its journey.

From Figure 15 we see that there is a clear indication that the earlier the cluster has its first GMC encounter, the higher $f_{\odot}$ it ends up with in the end. This is reasonable since the cluster will have a cold orbit in the beginning of its lifetime. Stars that escape early will therefore have a higher probability of ending up on a cold solar-like orbit. Additionally, clusters initially on dynamically hot orbits spend less time in the disc and hence have later first encounters.

The diffusion of the orbital elements of stars that escape cluster A during direct GMC encounters can be seen in Figure 8. This tells us that the nature of the diffusion is controlled by the galactic gravitational field and that it is only possible to link stars to a specific GMC encounter if the encounter has happened within a few galactic orbits. Stars that escape clusters that are located on high $|z|$ orbits stay in similar orbits as the cluster for a relatively long time since the influence of spiral arms is weaker and GMC encounters only happen close to the Galactic plane.

The average stellar escape rate between GMC encounters for clusters A, B and C ranges between 4.5 to 6.4 stars per Myr which is quite different from what Hurley et al. (2005) found in their $N$-body simulation of M67. Their average escape rate is $\sim 8$ stars per Myr for a cluster with 36000 initial stars with a binary fraction of 50 % and $r_{h}=3.9$ pc. Their cluster is placed at a Galactocentric distance of 8.0 kpc in a static galactic potential which makes our representation of M67’s galactic enviroment more realistic. The discrepancy between our clusters and those of Hurley et al. (2005) is dominated by our lack of primordal binaries. This makes the relaxation time for Hurley et al. (2005)’s cluster shorter, leading to a higher escape rate. It will also make the cluster dynamically older and is why in terms of remaining stars after 4 Gyr our clusters are somewhat similar with their cluster having 3520 stars left (870 singles and 1325 binaries). Our galactic enviroment and the GMC encounters for cluster A, B and C are in this way compensating for the lack of primordal binarity in our clusters in terms of ending up with the similar number of stars as Hurley et al. (2005). If we would want to have a more realistic dynamical version of M67 and add a binary fraction of 50 % we would need to add an additional number of $\sim 10^{5}$ stars to our initial clusters since the shorter relaxation time would increase the average stellar escape rate. A more massive birth cluster representing M67 could possibly explain the difference in the observed number of around 30 blue stragglers (Deng et al., 1999) versus the 20 found in the $N$-body simulations of M67 by Hurley et al. (2005). A more massive cluster would produce more blue stragglers and mass segregation in the cluster would cause the blue stragglers to sink to the cluster centre where they will be less affected by the galactic tidal forces. This would inadvertently increase the fraction of blue stragglers in the cluster as lower mass stars are tidally stripped from the cluster. We chose not to use binaries in our $N$-body simulations of M67 as a consequnce of computational issues with the binaries not working well with the tidal tensor of nbody6tt. In turn for using nbody6tt and no binaries we get a better representation of the galactic tidal field affecting the clusters, while only lowering the stellar escape rate which should not have significant effect on $f_{\odot}$.

7 Conclusions

In this paper we have simulated the evolution of M67-like clusters in a galactic Milky Way potential, followed the galactic orbits of escaping stars and compared them to the nature of the current solar orbit. We find that the simulated M67 candidates can be divided up into three groups; hot clusters, depleted clusters and scattered clusters. The hot clusters are born with $z$-velocities from the tails of the initial distribution. They therefore have an easier time reaching high $|z|$ without having to go through as many GMC encounters. Depleted clusters and scattered clusters all start out with a low $|v_{z,0}|$ and therefore need more heating by GMC encounters to reach the present day orbit of M67. Depleted clusters go through the largest number of GMC encounters, between 8 and 10, and as a result are depleted of stars to the point where they do not represent M67 at the present day.

All model clusters have stellar escapers that end up on solar-like orbits. Hot clusters have the lowest fraction of escapers that end up on solar-like orbits ($f_{\odot}=0.06$ %) and depleted clusters have the highest with $f_{\odot}=6.61$ %. Scattered clusters lie in between with an average $f_{\odot}$ of 2.1 %.

Hot clusters are similar to the clusters simulated by Pichardo et al. (2012) and we show that it is in fact possible to link the dynamics of these clusters to escapers with solar-like orbits by including GMCs to our galactic model. This dynamical link does not require the stars to escape the cluster with high escape speeds as seen in Figure 17, with the majority of stars having escape speeds less than 10  km s*-1*. These escape speeds are compatible with the survival of the wide planetary objects in the Solar System for all types of clusters. Stars above $\sim 10$  km s*-1*  are either black holes or neutron stars that have received kicks. The depleted clusters and scattered clusters correspond to the clusters that are investigated by Gustafsson et al. (2016) where the nature of scatterings between GMC and clusters was explored.

We have shown that the orbits of stellar escapers in a time-varying galactic potential diffuse from their original values over a few galactic orbits which means that only recent stellar escapers will be able to be dynamically linked to their host cluster.

By using $f_{\odot}$ from the depleted clusters and scattered clusters we make a rough estimate and calculate that clusters that are destroyed in the Galactic disc have a specific frequency of escapers that end up on solar-like orbits that is $\sim$ 2 times that of escapers from clusters that survive their journey.

Acknowledgements

We would like to thank Bengt Gustafsson, Melvyn Davies, Florent Renaud, Guido Moyano and Lennart Lindgren for valuable discussions and suggestions that helped form this paper. This work has been supported by grant 2017-04217 from the Swedish Research Council and the project grant 2014.0017 ‘IMPACT’ from the Knut and Alice Wallenberg Foundation. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc, which we can contribute thanks to grants from The Royal Physiographic Society of Lund.