© ESO, 2017

1. Introduction

The formation of the Oort Cloud is one of the important issues when trying to understand the origin and evolution of the Solar System. This has been the case ever since this structure was first recognised (Oort 1950), and resolving the issue still presents a very difficult task. It is natural to think of a “primordial” origin connected to the formation of the planets during the earliest stages of the Solar System more than 4.5 Gyr ago, as did Oort himself, and this has led to the classical picture (Duncan et al. 1987; Dones et al. 2004) of comets as icy planetesimals scattered through the gravity of the growing giant planets into orbits extending far enough to sometimes be decoupled from the planetary system by external agents (Galactic tide and passing stars).

A different scenario for the Oort Cloud formation was recently investigated by Brasser & Morbidelli (2013). They explored one of the consequences of the Nice Model (Levison et al. 2011, and references therein) for the long-term dynamical evolution of the giant planets. As a result of the rapid migration of Uranus and Neptune through the primordial trans-planetary disk into their current orbits, a scattered disk would be formed and, as an unavoidable by-product, also an Oort Cloud. This suggestion places the origin of the cloud at the time of the Late Heavy Bombardment (LHB) about 4 Gyr ago.

As long as this version of the Nice Model stands, there is strictly no need for the Oort Cloud to include any primordial component. However, even so, such a component is not ruled out. Moreover, the Nice Model may also accommodate a different scenario, where the planet migration happens very early. In this case, the Oort Cloud would definitely be primordial, and hence this option needs to be considered. An important issue then concerns the efficiency in the transfer of planetesimals into the Oort Cloud. There are two steps involved: first, the scattering of planetesimals into orbits that may be modified by external actions, and second, the decoupling that causes storage into the Oort Cloud.

It has been realised – ever since the work of Gaidos (1995) and Fernández (1997) – that the Oort Cloud storage is strongly dependent on whether one treats the new-born Sun as an isolated star or as a member of a dense stellar environment in a so-called birth cluster. The latter offers a more efficient way to decouple the objects by raising their perihelion distances due to the frequent occurrence of close and slow stellar encounters.

A numerical study of the formation of the Oort Cloud in a stellar cluster was performed by Fernández & Brunini (2000). In this work, the cluster was assumed to exist for a period of 100 Myr with a stepwise decreasing number density of stars from the assumed initial value down to zero. In addition, the tidal effect of the placental molecular cloud gas was included, typically only for the first 10 Myr. The scattering and decoupling of comets was simulated with the main result being that an Oort Cloud inner core was formed quite rapidly with Jupiter and Saturn as the main scattering agents.

The discovery in 2003 of (90377) Sedna, whose perihelion distance of 76 au is well beyond the orbits of the giant planets and thus can only be explained by the influence of external actions, spurred an interest in very dense stellar environments for the birth of the Solar System. Embedded clusters were recognised to be a common birth place for solar-type stars (Lada & Lada 2003). Brasser et al. (2006) found a good match to the Sedna orbit for an inner Oort Cloud obtained in a model, where the Sun was born in such a cluster with a very high mass density. However, in a follow-up paper, Brasser et al. (2007) investigated the effects of gas drag from the solar nebula on the planetesimal scattering and found what they referred to as size sorting. Only very large objects would evolve in the way described by Brasser et al. (2006), while orbits of comet-sized objects (radii ~1 km) would be circularised beyond the planetary orbits, opposing their scattering.

In a paper by Kaib & Quinn (2008), the formation of the Oort Cloud in a stellar cluster environment was again considered with similar assumptions for the cluster lifetime and number density of stars. However, the evolution thus computed during the first 100 Myr was supplemented by an additional 4.4 Gyr in a model of the galactic disk tide plus field star encounters. The effect of the resulting present Oort Cloud being dominated by a tight, inner core was confirmed as well as the possibility of obtaining a Sedna-type population due to the random effects of the closest stellar encounters.

An often cited model for the formation of a primordial Oort Cloud was proposed by Levison et al. (2010). This model assumes a birth cluster with few stars (30 <N< 300) and thus with a short lifetime. In this case, the stars are found to fly apart when the gas component of the system is purged due to external influences. The Oort Cloud is formed by scattered disk comets from different stars of the cluster and becomes enriched during the cluster break-up. However, the above-mentioned problem of bringing kilometre-sized objects into extended scattered disks was not addressed.

The current models of Oort Cloud formation may thus be summarised as, on the one hand, models for a primordial cloud, assumed to be formed in a dense but more or less short-lived stellar environment, and on the other hand, a delayed formation model associated with the LHB, where the Solar System is assumed to have left its birth cluster at an earlier stage. All these models assume a rather short dissolution lifetime for the birth cluster, which may or may not reflect the truth.

The size of the Sun’s birth cluster was discussed in a review by Adams (2010). His analysis of a range of constraints led to a broad probability distribution for the number N of stars, peaking at N ≃ 2500. An important argument against overly small values of N was that it would be too unlikely for the birth cluster to produce a supernova with a progenitor mass of 25 ℳ_⊙ or more, as seems necessary to explain the amounts of short-lived radio isotopes that meteorite evidence show to have been present in the solar nebula (e.g. Williams & Gaidos 2007). However, in itself, this factor gives rather weak constraints – only clusters with N < 50 would be excluded, since the random likelihood of the supernova would then be less than 5%. The main factor opposing overly large N values was an overly small chance for the regularity of the giant planet orbits to survive in the presence of the close stellar encounters then implied (Adams & Laughlin 2001; see also Malmberg et al. 2007).

More recently, Gounelle & Meynet (2012) proposed a model for the origin of the short-lived radionuclides involving two generations of stars, formed in the same giant molecular cloud and preceding the formation of the Sun. This model implies a Solar birth environment rich in stars. Several thousand stars were estimated to have been born before the Sun, thus providing the source for ²⁶Al and ⁶⁰Fe traced in chondritic meteorites. The dynamical fate of the Sun was not addressed, but it seems possible that the Sun stayed gravitationally bound to the initial complex of stars and gas, thus becoming a member of a massive stellar cluster.

Thus, the Sun’s birth cluster may have been rich in stars, containing thousands or more to begin with. Moreover, if we leave aside the specific constraints posed by solar system evidence and consider the statistics of observed embedded clusters, we find that the number of clusters formed today falls off with the number of member stars in such a way that approximately equal numbers of stars form in clusters with 10², 10³, and 10⁴ members (Lada & Lada 2003). While several papers have dealt with Oort Cloud formation in a cluster with ~10² members, which we shall call a low-mass (LM) cluster, the other two classes of clusters – intermediate-mass (IM) and high-mass (HM) clusters – have not yet been considered. The cases treated by Fernández & Brunini (2000) and Kaib & Quinn (2008), where the whole cluster dissolves within 100 Myr, would fall into our LM category.

Dynamical models of stellar cluster evolution show that clusters with more than 1000 initial members typically survive for several hundred Myr or more (Lamers & Gieles 2008). Thus, their lifetimes may even exceed the interval from the formation of the earliest solar nebula condensates (meteoritic calcium-aluminium rich inclusions or CAI) until the triggering of the LHB (Morbidelli et al. 2012). This calls for a reevaluation of the Oort Cloud formation models in the framework of such a long-lived birth cluster. The large number of stars may help in the formation of a primordial cloud, but it is also a threat to the stability of the cloud due to the possibly disruptive effect of subsequent encounters. In the present paper we consider the fate of a primordial Oort Cloud in a dense stellar environment that lasts at least until the LHB. The main question is if such a cloud would survive or not.

We use two assumptions for the birth cluster, considering two values for the initial number of stars (N₀). For the IM cluster we take N₀ = 2000, and for the HM cluster we choose Messier 67 as a template. In this case, N₀ = 36 000 whereof half the systems are binary (Hurley et al. 2005). This choice is arbitrary, and we do not mean to suggest M 67 to be the Sun’s birth cluster – for discussions of this issue, we refer to Pichardo et al. (2012) and Gustafsson et al. (2016). It is, however, a convenient case for the upper end of the mass range due to the availability of detailed evolutionary modelling based on a very good observational record.

We have devised a modelling technique that allows us to trace the motions of thousands of test objects representing Oort Cloud comets on heliocentric orbits for a time span of several hundred Myr, in the presence of a static cluster potential plus a selection of randomly created stellar encounters. These are meant to include some of those that impart the largest impulses to the Sun. The CPU time consumption is moderate enough to allow running thousands of such integrations for each model of the birth cluster and thus obtaining results that are robust against statistical uncertainties.

In Sect. 2 we first present our cluster model and describe its two versions (IM and HM) in some detail. We then describe our treatment of stellar encounters and our derivation of the encounter frequency and velocity distribution. Next, we present our simulation set-up. The results are given in Sect. 3. Notably, we find that the survival probability of primordial Oort Cloud comets is extremely low regardless of the size of the birth cluster, within the cluster mass range explored. Conclusions drawn from this result and a discussion are given in Sect. 4. In two appendices we describe the calculation of the cluster model structure and the implementation of the stellar encounters, respectively.

2. Methods

We integrate the orbit of the Sun together with thousands of test particles representing Oort Cloud comets, in a static cluster potential over the course of 400 Myr. The comets are introduced in random heliocentric orbits with semi-major axes of a_o = 5000, 10 000, and 20 000 au. The cluster potential is computed for an intermediate-mass (IM) cluster based on a template with an initial number of N₀ = 2000 stars, and a high-mass (HM) cluster whose template initially had N₀ = 36 000 stars, as described in Sect. 2.1. We represent the HM cluster by one static model while we use a sequence of static models for the IM cluster, as this evolves rapidly with time and in fact nearly dissolves during the simulated time period.

During each simulation, we select roughly 20 stellar encounters to occur at random times, by means of an impact approximation described in Sect. 2.2. We add these interloping stars one at a time to the orbit integrations as described in Sect. 2.3. The gravitational influence on the comets is thus that of the smooth cluster potential, the Sun, and at times, one additional star. The numerical setup and the selection of comets is further described in Sect. 2.4.

2.1. General properties of the cluster model

We impose the following constraints on our cluster model.

• The stellar distribution function corresponds to a spatiallyisotropic, relaxed state.

• It does not directly account for mass segregation and considers only single stars.

• It is stepwise or fully constant with time during the interval we simulate.

These constraints are to some extent mutually incompatible, since mass segregation and binary formation are necessary consequences of the same dynamics that causes relaxation. In addition, for the relatively young system that we simulate, relaxation should still be ongoing. Thus, the system has to evolve with time in contrast to our third constraint. The reason why we stick to the above concepts is that they allow us to develop a synthetic model of stellar encounters that greatly facilitates the Oort Cloud simulations and avoids the use of time-consuming N-body simulations.

Clearly, a most useful template in order to achieve a cluster model with the above properties is the King model (King 1966). We apply this concept according to a prescription that we describe in Appendix A.

To fit our model parameters, we start from an initial guess (see Appendix A), and improve these values iteratively, until they yield the desired solution. This means that the cluster has the required mass, and the density function drops to zero at a distance R_lim less than the tidal radius. Other fitted parameters include the half-mass radius r_h. Finally, we may also compute the sky projected surface density I(R) as prescribed by King (1966, Eqs. (23)–(27)). Using this quantity, we calculate the core radius R_c from its definition $\mathit{I}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{I}\mathrm{\left(}\mathrm{0}\mathrm{\right)}$. In our case, the surface density refers to mass rather than brightness. This core radius may also be used as a fitting parameter.

The last property to be defined is the age of the cluster. We want to simulate the cluster effects on the Oort Cloud for a time interval, extending from the formation of a primordial cloud until the LHB. However, in this work we do not consider the effects of the cluster on the very formation of the cloud. In reality, the formation process is expected to extend over several hundred Myr (Kaib & Quinn 2008; Brasser & Morbidelli 2013), but we replace this by a step function: we consider a first time interval of 100 Myr starting at the formation of the solar system (t = 0), during which there is no Oort Cloud at all. Then, at t = 100 Myr, we introduce the comets in their above-described, initial orbits. The cluster age hence extends from 100 Myr to 500 Myr, when we assume the LHB to occur.

2.1.1. The high-mass (HM) cluster

As indicated above, the structural parameters of this cluster are based on the properties of a young M 67 as simulated by Hurley et al. (2005) and available in their online tabulations¹. These authors performed a set of N-body simulations, describing the complete evolution of M 67 after formation of all the stars and escape of any residual gas. The initial parameters of their evolutionary model – the galactocentric distance and orbit, total mass, and binary populations – were estimated from the current luminous mass (i.e. the total mass of nuclear-burning stars, estimated at ℳ_L ~ 1000 ℳ_⊙), age (approximately 4.0 Gyr), and binary and blue straggler populations of the cluster. As these populations are sensitive to the dynamical evolution of the cluster, the model is reasonably well constrained.

In Fig. 1 we present the evolution of both our template clusters. The HM cluster, based on data from Hurley et al. (2005), is seen to decrease slowly in the number of stars during the considered time interval (from 100 to 500 Myr). We illustrate the total number of stars, which often occur as binary components in the real cluster M 67 – a discussion of binarity is given in Sect. 2.1.3. The inner part of the cluster undergoes a slight expansion with the half-mass radius increasing from an initial value of 4.0 pc to 5.5 pc.

Fig. 1 The number of stars (solid, black lines) in our high-mass and intermediate-mass cluster models, and their half-mass radii (dashed, red lines), versus time from cluster formation. A logarithmic scale is used for the number of stars and a linear scale for the radius. The time has been rescaled for the IM cluster to take the expected influence of GMCs into consideration (see text).

The initial distributions of mass and energy used by Hurley et al. (2005) were based upon a Plummer model (Plummer 1911). On the half-mass relaxation timescale T_rh,0 ≈ 290 Myr, this distribution evolves into something resembling a King profile. The relaxation causes the increase in the half-mass radius, which reflects the evolution of the density profile.

We choose a static setup for our model cluster representing the mean state of M 67 during the age span 100–500 Myr in the simulations of Hurley et al. (2005). At this early stage in the cluster’s history, the state can be described as semi-relaxed, so that the global structure is neither that of the initial Plummer distribution nor that of the eventual King model. However, the relaxation proceeds more rapidly in the inner parts of the cluster, where the crossing time is smaller. Therefore, the structure in this part of the cluster is less sensitive to the initial conditions of their model. By choosing a King model, we are thus able to reproduce the inner parts of the simulated cluster reasonably well, while the outer parts are less well described.

The parameters to which we fit our King model are the mass ℳ_cl, half-mass radius r_h, core radius R_c, limiting radius R_lim, and the mean density within the half-mass and core radii, ⟨ ρ_h ⟩ and ⟨ ρ_c ⟩. The input parameters and resulting properties of the computed model are given in Table A.1.

2.1.2. The intermediate-mass (IM) cluster

While N-body simulations of the HM cluster indicated a rather slow evolution during the first 2–3 Gyr, the IM cluster has a short relaxation time and thus rapidly evolves not only in total mass but also in terms of structure. The half-mass relaxation time scales with the cluster mass and half-mass radius like ${\mathit{t}}_{\mathrm{rh}}\mathrm{\propto }{\mathrm{ℳ}}_{\mathrm{cl}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{r}}_{\mathrm{h}}^{\mathrm{3}\mathit{/}\mathrm{2}}$ (Spitzer & Hart 1971). For a population of N₀ = 2000 stars with r_h = 1 pc, we find t_rh ≈ 20 Myr, indicating that such a cluster should relax very rapidly.

In Fig. 1 we illustrate the simulated evolution of an IM cluster, computed using the emacss code (Alexander & Gieles 2012; Gieles et al. 2014; Alexander et al. 2014), with its lifetime rescaled according to the expected influence of encounters with GMCs (Gieles et al. 2006). We find that the relaxation time varies between 20 and 40 Myr during the first 500 Myr of the simulation. The parabolic evolution of the half-mass radius after the initial relaxation is due to the combined effect of mass loss from the outer regions, and the collapse and bounce of the core, where energy transfer from the core inflates the outer regions of the cluster.

Due to the rapidly evolving nature of this cluster, we represent it by a sequence of static models, each corresponding to the average properties during a span of 25 Myr. The input parameters and resulting properties of the models are listed in Table A.1. The density distributions of the HM cluster model as well as every fourth model representing the IM cluster, are shown in Fig. 2. Both structures exhibit an inner plateau out to the core radius (1.9 pc in the HM cluster; 0.05–0.1 pc in the IM cluster), whereafter the density roughly follows a power-law, decreasing to values lower than that of the local Galactic disk.

Fig. 2 Density distribution in the template clusters. The dashed black curve indicates the distribution of the high-mass cluster, while the solid coloured curves represent time steps in the intermediate-mass model, 100 Myr apart, with the average time of each step shown in the legend. The central mass density (left end of the plot) decreases monotonically with age of the model. The horizontal dash-dotted line indicates the Galactic mid-plane density for comparison.

2.1.3. Cluster stars

In addition to the structure of the cluster, we need to describe its stellar content. To translate a distribution of mass into the corresponding distribution of individual stars, we adopt a stellar IMF from the generating function of Kroupa et al. (1993), $\mathrm{ℳ}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{=}{\mathrm{ℳ}}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{0.19}{\mathit{\xi }}^{\mathrm{1.55}}\mathrm{+}\mathrm{0.050}{\mathit{\xi }}^{\mathrm{0.6}}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\xi }{\mathrm{)}}^{\mathrm{0.58}}}\mathit{,}$(1)where ξ ∈ [0,1] is a random number from a uniform distribution, and ℳ₀ the lower mass limit. We set ℳ₀ = 0.1 ℳ_⊙ to account for the mass segregation-driven preferential loss of the lowest-mass objects from the cluster.

We evolve the stellar population to an age of 300 Myr, representing the mean age of the cluster, using the rapid stellar evolution code SSE (Hurley et al. 2000). The resulting mean stellar mass is ⟨ ℳ ⟩ ≈ 0.41 ℳ_⊙, with an effective upper mass limit near 3.5 ℳ_⊙ (turnoff mass ~3.3 ℳ_⊙). For completeness, we retain stellar remnants in the form of white dwarfs in the mass distribution, but remove neutron stars and stellar mass black holes as these are expected to usually be given kicks at formation much greater than the cluster escape velocity (see, e.g. Pfahl et al. 2002).

Fig. 3 Distribution of relative velocities in our encounter model for the stellar cluster. The five curves in the left panel represent, from top to bottom, the maximum allowed relative velocity (solid, black), its 90th percentile (dotted, red), mean (solid, blue) and 10th percentile values (dotted, black). Shown as vertical (black) lines are the core (dotted) and half-mass (dashed) radii. Left panel: relative velocity distributions for the high-mass cluster. Right panel: relative velocity distributions for the intermediate-mass cluster at two time steps. The first step represents the time interval 100–125 Myr. The second time step represents the time interval 300–325 Myr, and is shown only by a thick dashed black line representing the maximum relative velocity, a thick dashed blue line representing the mean, and a vertical thick dashed line representing the half-mass radius at this time. The core radius is not shown at this time, but is available for all time steps in Table A.1.

We do not directly invoke mass segregation in our cluster model. Thus, the probability distribution of stellar mass is the same at any distance from the centre. In this sense, our model fails to include a real phenomenon, which would be expected to occur in either cluster a few hundred Myr after its birth; a systematic tendency for the high-mass stars to concentrate in the cluster core and, hence, to be underrepresented in other parts of the cluster.

We also neglect the existence of binaries, even though these are common in M 67 (Richer et al. 1998) and also in low-mass clusters (Giersz & Heggie 1997), since they would complicate the treatment of close encounters. If hard binaries were included, these would effectively increase the mean mass of the Oort Cloud perturbers while for a given cluster mass decreasing the overall number density of those perturbers. We note that Fouchard et al. (2011) found the long-term dynamics of the Oort Cloud to be influenced by massive stars in the Galactic disk to a much larger extent than their low encounter frequency would suggest – the reason being their higher chance of producing global perturbing effects on the cloud. Hence, we may be underestimating the stellar encounter effects, but this is in line with our strategy to seek a conservative estimate.

2.2. Stellar encounter flux

Since we do not trace the motions of individual cluster stars, we generate encounters synthetically using a statistical encounter flux derived from the cluster model.

To derive this encounter flux, we first need the distribution of relative velocities of the stars. At a given distance from the cluster centre, we calculate the cumulative distribution of kinetic energies from Eq. (A.1) and produce a generating function. To represent individual encounters, we draw independently two random values from this generating function for the kinetic energy (per unit mass) of two encountering stars.

By assuming a flat distribution of angular momenta (as inherent in a King model), absorbing, at most, all of the available kinetic energy, we determine the tangential and radial velocity components of each star. Each radial velocity can be positive or negative with equal probabilities. We assume an isotropic encounter distribution, wherein the motions of the two stars are uncorrelated. The distribution of these encounter velocities is illustrated in Fig. 3, including the arithmetic mean ⟨ v_rel ⟩ which in the figure is denoted v_mean and shown by blue curves. We note that when we use this to find the Sun’s encounter flux, we do not account for the constraints imposed by the Sun’s particular orbit. In our approximation, the Sun is treated as an average star at any particular distance from the cluster centre.

Fig. 4 Encounter frequency versus distance from the cluster centre, shown for the high-mass cluster (dashed line) and selected time steps for the low-mass cluster (coloured solid lines). The encounter flux at the centre of the cluster decreases with increasing age for the intermediate-mass cluster – colours are the same as in Fig. 2. The horizontal dash-dotted line illustrates for comparison the current encounter frequency of the Sun with Galactic field stars (Rickman et al. 2008).

Using the mean relative velocities, we show in Fig. 4 the encounter frequency versus distance from the cluster centre: f(r) = n ⟨ v_rel ⟩. For an assumed impact parameter b, this encounter frequency produces an encounter rate Γ ≡ σn ⟨ v_rel ⟩, where σ ≡ πb² is the impact cross-section. The expected number of encounters occurring over time τ with a constant value of the encounter rate is thus simply ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\tau }}\mathrm{\Gamma d}\mathit{t}\mathrm{=}\mathit{\tau }\mathrm{\Gamma }$. Taking Γ = τ^-1 for a given time τ, we may thus estimate the expected minimum impact parameter b_min of the encounters experienced during a timespan of length τ. Using τ = 400 Myr, in the HM cluster, we find a value of b_min ≈ 210 au at the very centre, while near the half-mass radius we find b_min ≈ 1200 au.

For the IM cluster, we integrate ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\tau }}\mathrm{\Gamma }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}\mathit{t}$, taking into account the time dependence of the encounter flux, by adding the contributions from the different 25 Myr intervals. At the cluster centre, we find b_min ≃ 35 au for τ = 400 Myr. Clearly, and not surprisingly, the Oort Cloud and even the planetary system would not survive such an extended stay in such a dense environment. If we take into account that the Sun is found to leave this kind of cluster before 100 Myr as a median (see Sect. 3.4), for this value of τ we still have b_min ≃ 40 au, which is very destructive. However, the solar orbit makes the Sun spend only a minimal amount of time in the immediate vicinity of the centre, if any at all. For a more realistic estimate, we also perform the estimate for a distance of 1.1 pc from the centre, roughly representing an average half-mass radius. In this case, the two values of b_min increase to approximately 1450 au and 2000 au, respectively. These are larger than the one found for the HM cluster.

The expected minimum impact parameter resembles the size of the Oort Cloud, b_min ≈ 10 000 au at r ≈ 14 pc in the HM cluster, or at r ≈ 3.5 pc in the IM cluster. We note, however, that at these distances in either cluster, the cluster mass density is in fact comparable to that in the Galactic disk (~0.1 ℳ_⊙ pc^-3, Holmberg & Flynn 2000). In the disk, the mean encounter velocity is greater by an order of magnitude, which increases the flux of stellar encounters by the same amount to 10^-5 pc^-2 Myr^-1, but also reduces the efficiency of momentum transfer. The influence of field star passages on the Oort Cloud has been investigated elsewhere (Rickman et al. 2008), and the erosion of the Oort Cloud was then found not to be dramatic over intervals like the one considered here. Thus we neglect these passages as well as their associated Galactic disk tide effect.

2.3. Stellar encounter selection

We simulate stellar encounters by introducing a star at a distance d_start from the Sun, and integrate the orbits of the two stars under their mutual gravitational influence in the cluster gravitational potential, until the mutual distance is d_end. In the simulation of the HM cluster, we found that setting d_start = d_end = 1 pc results in a typical encounter duration less than 1 Myr, with a tail in the distribution of encounter durations extending toward 16 and 34 Myr representing fewer than five and one percent of encounters, respectively. Such extended encounters are similar in duration to the period of a Kepler orbit around the cluster centre near the half-mass radius in this cluster, $\mathit{P}\mathrm{~}{\mathit{r}}_{\mathrm{h}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{ℳ}}_{\mathrm{cl}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{~}\mathrm{9}\mathrm{Myr}$.

In the simulation of the IM cluster, a distance d_start of 1 pc would be inappropriate as this is in fact similar to the half-mass radius of the cluster. Thus, most of the interaction of the two stars occurs at distances much larger than those between the Sun and many other cluster stars, which makes the calculation somewhat irrelevant. In addition, as we shall explain below, we use a filter when selecting encounters to avoid overlaps of consecutive events. With the smaller relative velocities of the stars in the IM cluster, this would lead to the blocking of too many low-velocity encounters due to their large durations. To reduce this bias, we set a distance d_start = 0.5 pc where the encountered star is introduced, and a slightly smaller distance d_end = 0.45 pc where it is removed. This results in a similar distribution of encounter durations as in the case of the HM cluster.

We have checked that the particular choice of distances does not significantly influence the results by repeating our IM cluster calculations with d_start = 0.25 pc and d_end = 0.22 pc, with no significant effect on either the survivability of comets or the evolution of the solar orbit. Similar tests on the HM cluster indicate that neither qualitative nor quantitative results depend strongly on the precise choice of the interaction distance.

The typical encounter in either simulation does not overlap subsequent pericentre passages, which is important as our setup does not allow for simultaneous encounters. We enforce this by a veto blocking the selection process while an encounter is ongoing. To avoid having the encounter scheme block too many subsequent encounters (typically occurring near pericentre), we aim for producing one encounter per 20 Myr, i.e. Γ = (20 Myr)^-1. For a simulation duration of τ = 400 Myr, this corresponds to typically 20 encounters per simulation, where the veto typically blocks one expected encounter per simulation.

The number of stellar encounters expected to exert significant influence on the outer Oort Cloud is, however, larger than just 20 per 400 Myr, and we aim to simulate those encounters that are most important. Rather than using the distance as criterion, we adopt a strength parameter S = ℳ/(v_relb), approximating the impulse transferred to the Sun by a stellar encounter in the classical impulse approximation (see Rickman 1976; Fouchard et al. 2011). All three defining parameters must then be known for S to be computed, which in turn requires the encounter geometry to be determined. As detailed in Appendix B, we pick the encounter parameters at random and generate a list of encounters likely to occur during a time interval of a given length, and compute S for each. Finally, we select the encounter with the largest value of S.

When the final selection is made, the integration proceeds until the encounter is finished, as detailed in Sect. 2.4. Then a new time interval of the given length is considered, using the modified solar orbit, and the next stellar encounter is selected. On average, the interval between consecutive encounters will be close to 20 Myr, and hence a total of about 20 encounters will be treated during the entire 400 Myr simulation.

2.4. Numerical simulations

We use a hybrid numerical setup. We combine the tidal effects of the cluster as a whole, represented by the smooth gravitational potential computed according to Sect. 2.1, with the influence of individual stars during orbit-integrated close encounters as described in Sect. 2.2. Comet motions are thus integrated under the gravitational influence of the cluster, the Sun, and an interloping star when applicable. As mentioned above, we neglect the influence of Galactic tides and interloping Galactic field stars.

The Oort Cloud is represented by a sample of 3000 comets. This consists of three ensembles of 1000 comets each, representing different parts of an initial Oort Cloud, with original semi-major axes of a_o = 5000, 10 000, and 20 000 au. We shall here refer to the respective (initial) ensembles as the (initially) inner, intermediate and outer comet cloud – not to be confused with the nomenclature of actual present-day Oort Cloud populations. The other orbital elements are drawn randomly from identical distributions for all three ensembles. These distributions are uniform for all but the eccentricity, which is distributed as f(e) ∝ e in the range e ∈ [0,1] to represent a thermalised state. Since our model does not include the planets, the initial cloud is modelled without any loss cone (Oort 1950; Hills 1981), thus allowing e → 1. Orbital inclinations are given an isotropic distribution, uniform in cosi ∈ [−1,1].

Fig. 5 Solar orbit during the first 80 Myr in a simulation of the high-mass cluster, with grey circles indicating the core and half-mass radii. The solar orbit is illustrated by a line that switches colour (from violet to blue, green and yellow) every time a star is encountered. The starting points of the interloping stars are numbered and indicated by a star symbol, while end points are indicated by x. The starting point of the solar orbit is indicated by a circle of radius 20 000 au. The solar orbit is by construction initially confined to the x-z plane, and departs significantly from that plane only after a close interaction with star number 3.

The simulations are realised in the barycentric frame of the cluster, with the Sun initially positioned at the half-mass radius in the case of the HM cluster. In the IM cluster, we instead placed the Sun at a random distance from the cluster centre, distributed according to the cluster density profile. The initial velocity (and orbit) of the Sun are determined randomly from the energy and angular momentum distributions, and the time and configuration of the first (or next) encounter are determined as was described in Sect. 2.3.

The orbit of the Sun is integrated together with the comets until the initialization of the following encounter, where the star is introduced at a distance of 1 pc or 0.5 pc from the Sun for the HM or IM cluster, respectively. The integration then includes the mutual gravitational influence of the Sun and the star on each other, and on the comets, until the distance between the stars is again 1 pc or 0.45 pc for the HM or IM cluster, respectively. All the objects are also subject to the smooth gravity field of the cluster. An example of the solar orbit during the first 80 Myr of one of the simulations of the HM cluster is given in Fig. 5.

Comets that move beyond a distance of 1 pc from both the Sun and the star are logged and then discarded from the computations. In addition, following the results of Fouchard et al. (2013), we model a planetary loss cone such that comets reaching within 5 au of the Sun are considered ejected by Jupiter. For simplicity, the same condition is applied when the comet is close to the interloping star, assuming it is orbited by a similar giant planet.

Numerical integrations are performed with the 15th-order RADAU integrator of Everhart (1985). With the limited time span of our integrations compared to the time step used, this integrator is perfectly suitable. It is not energy preserving (symplectic), but for the problem at hand, small departures from energy conservation are not an issue, and the integrator has been found to perform very well in closure tests (Carusi et al. 1985).

3. Results

We perform a series of 1000 Monte Carlo simulations for each cluster, differing from one another in terms of the initial solar orbit and cometary orbits. The total number of comets treated is thus 3 × 10⁶, for each cluster. With such statistical sampling, we shall present not only the typical behaviours, but also the tail of the probability distribution, that is, the rare outcomes. Obviously, in the latter cases, our results are not statistically accurate but rather indicative.

3.1. Comet survival probability

Fig. 6 Distribution of comet survival probability for the three different bins of initial semimajor axis, ao = 5000 au (red), 10 000 au (blue), and 20 000 au (orange). Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster.

We illustrate in Fig. 6 the statistical distribution of comet survival probability in these simulations. The term “survival” refers to those comets that do not experience ejection either by moving to distances exceeding 1 pc or by intruding into the loss cone – hence they stay in the Oort Cloud until the end of the integration. Since it is natural to imagine that stellar encounters cause comets to leave the cloud directly, while the cluster tide may perturb their orbits into the loss cone, we have also performed two extra sets of 100 simulations each, where we artificially turned off the cluster tide and the stellar encounters, respectively. By comparing these results with those of the full model, we hope to learn which is the dominant mechanism causing the loss of comets from the Oort Cloud.

Results for the HM cluster are presented in the left panel of Fig. 6. In this case, the median survival probability for comets in the initial sub-population representing the primordial inner comet cloud, a_o = 5000 au, is just 0.2%, representing a few comets per simulation. For comets representing the primordial intermediate and outer comet clouds, a_o = 10 000 au and 20 000 au, the median survival probability is less than 0.1%, that is, fewer than one comet per simulation. Hence, typically the entire primordial outer comet cloud is lost. In only 1.3% of simulations do more than 1% of the primordial outer cloud members survive. For the primordial inner and intermediate comet cloud, more than 5% of comets survive in 9% and 2% of the simulations, respectively.

Comparing now with the results for the IM cluster in the right panel of Fig. 6, we see an important difference. The curves representing the fall-off of the percentage of simulations with increasing percentage of survivors are much flatter for the IM cluster over almost the whole range. Consequently, the median survival probability for the primordial inner comet cloud is only 0.1%, that is, even lower than for the HM cluster, while the percentage of simulations with a much larger number of survivors is considerably higher for the IM cluster. For instance, the fraction of simulations with more than 1% survivors in the primordial outer cloud is 9.8% for the IM cluster, compared to 1.3% for the HM cluster. On the other hand, the escape times of comets from the Sun in the two clusters are much shorter in the IM than in the HM case. For the inner, intermediate and outer clouds, the median escape times are 12, 5 and 4 Myr, respectively, in the IM cluster while in the HM cluster these are 51, 47 and 20 Myr.

The reason for these differences has to do with the typical fate of the Sun in the two clusters. As we shall see in Sect. 3.4, after 400 Myr the Sun typically remains in the HM cluster and escapes from the IM cluster. Specifically, the remaining percentage is 94% in the HM case and the escaping percentage is 95% in the IM case. This means that the behaviours exhibited in the two panels of Fig. 6 may actually carry as much information about whether the Sun remains in the cluster, as the difference between HM and IM clusters.

Although the statistics is rather poor for the less common situations, we can still compare the fate of Oort Cloud comets in all four cases, that is, HM versus IM clusters and remaining versus escaping Sun. We have thus found that for the IM cluster there is not much difference of comet survival statistics, whether the Sun stays in the cluster or it escapes. For the HM cluster the fate of the comets is more sensitive to the fate of the Sun. When the Sun escapes from this cluster, the comet survival statistics is intermediate between the left and right panels of Fig. 6.

Our results are consistent with the following picture. As seen in Fig. 4, the central part of the IM cluster starts out with a very high encounter frequency, but this falls off rapidly with time due to cluster evaporation. The trend is similar for the strength of the cluster tide. Thus, during the early phases the solar system runs a high risk of being stripped of its entire Oort Cloud due to very close stellar encounters, if the solar orbit penetrates close to the centre, but our simulation only covers a few of the strongest encounters expected. As we shall see below, the cluster tide also plays a role in this context. We may therefore expect to see a majority of disastrous cases with no or very few comets surviving and at the same time another category of cases, where the most perilous encounter was weaker and left a significant part of the Oort Cloud bound to the Sun. This situation quickly became fossilized, as the cluster started to dissolve and the Sun migrated outward before finally escaping.

Consequently, the time of escape might not matter very much for the survival statistics of the Oort Cloud, and even in case the Sun remains in the IM cluster for the whole interval considered, relatively little further damage to the cloud may be the rule. The case of the HM cluster is different, because its central region remains perilous for the full length of the integrations. Therefore, in case no disastrous encounter occurs during the early stage, the remaining part of the Oort Cloud will in general be subject to further damage due to close encounters. The chance for a significant fraction of surviving comets is relatively small. However, there are of course situations, where the Sun undergoes efficient outward migration, and the survival rate is higher. This will be the case in particular for simulations where the Sun escapes from the cluster.

As shown in Table 1, the comparison simulations where we artificially turned off either the cluster tide or the stellar encounters, indicate that the survival probability for the HM cluster is hardly affected at all by the presence or absence of the cluster tide in the dynamical model, while the removal of the stellar encounters drastically increases the chances of survival for all the primordial cloud populations. It is thus clear that the losses of comets are mainly due to the stellar encounters. However, for the outer cloud we also note that the median survival probability is very low even without stellar encounters. Hence, in this part of the cloud – and only in this part – the comets are destabilized also by the cluster tide.

We also see an indication in Table 2 that, for the IM cluster, both the cluster tides and the stellar encounters matter for the loss of comets from the Oort Cloud – the tides even more than the encounters. Here we have to take note of the fact that the scenario in the model without stellar encounters is very different from that of the other two models. In the full dynamical model the Sun escapes from the cluster in a large majority of cases, as seen above, and the same is true for the model without the cluster tide. However, when there are no stellar encounters, the solar orbit remains practically unchanged – hence, the Sun never escapes but stays close to its initial orbit during the whole integration. The initial orbits penetrate close to the cluster centre and the comets are hence sensitive to the radial tide. This explains the extensive losses of comets from all parts of the Oort Cloud.

3.2. Comet end states

In Tables 1 and 2 we present statistics regarding the three end states, for both clusters: the direct departures leading to unlinking of the comets, the entries into the loss cone leading eventually to hyperbolic ejection by Jupiter, and the survivals until the end of the integration. These are shown for the full dynamical model as well as the two comparison models, and for each set of primordial comet orbits. In each case, the listed percentages refer to the median of the simulations.

The most striking feature concerning the HM cluster (Table 1) is the predominance of direct departures in the full model as well as the model without tides. The model without encounters is different as regards the inner and intermediate comet populations, where instead of the predominance of departures we find large fractions of survival or loss cone intrusion. Once more we see that the outer population is very vulnerable to direct departures even without encounters. Apparently, the cluster tide causes a strong instability of the outer cloud orbits but much less so for the orbits of the other parts.

The loss cone entries practically only appear in the presence of the cluster tide, showing that stellar encounters very rarely cause such an evolution. On the contrary, stellar encounters are seen to interfere with the tidal evolution of the perihelion distance, preventing the loss cone entries from the inner and intermediate cloud that would otherwise occur. By plotting the perihelion distance versus time in the tide-only model, we have verified that the loss cone entries are caused by a secular oscillation of orbital angular momentum driven by the cluster tide. Even though such a pulsation is also a feature of the Galactic disk tide currently experienced by Oort Cloud comets (Heisler & Tremaine 1986), the dynamics inside the cluster is basically different. The cluster tide is radial and non-conservative. The amount of the energy exchange is shown by the preference for tidally caused departures of Oort Cloud comets belonging to the outer population.

Comparing with the results for the IM cluster (Table 2), the main difference appears for the model without stellar encounters. As noted above, the IM model is special in that the Sun remains more or less locked to its initial orbit for the full length of the integration. As discussed in Sect. 3.4, these orbits tend to have pericentre distances less than 0.5 pc. That this exposes the Oort Cloud comets to a very strong cluster tide can be realised from Fig. 2, because the strength of the cluster tide is proportional to the mass density in the homogeneous, central part. This density is seen to be very high at all times in the IM cluster, and the region of homogeneity extends to r ≃ 0.1 pc.

For this reason we see a very high fraction of unlinked comets in all parts of the Oort Cloud. Of course, the model without stellar encounters is not meant to be realistic. In reality the stellar encounters would rapidly change the solar orbit, thereby in general decreasing the influence of the cluster tide. What this model shows is that in an IM cluster the central region is very dangerous for Oort Cloud comets due to the cluster tide, and in case the solar orbit remains with a small pericentre distance for too long, the cloud is likely to be stripped away.

Fig. 7 Time evolution of the combined comet population of all simulations with the full dynamical model. The fractional populations are counted with respect to the total initial population. Results were sampled at 10 Myr timesteps, using 100 uniform logarithmic bins for a ∈ [500,200 000] au. The locations of the three initial sub-populations are indicated by vertical dashed lines. Here, snapshots are shown at times T = 10, 50, 200, and 400 Myr, from top to bottom, using black, red, blue and green colours, respectively, where the final time step is also shaded in grey. Comets beyond the plot limits have been summed at the edges. Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster. The three remaining peaks in the final distribution reflect the fact that the vast majority of the stars in this cluster model have by then left the cluster with fossilised structures of their cometary clouds.

3.3. Oort Cloud evolution

The time evolution of the comet cloud is illustrated in Fig. 7. We will first discuss the left panel, showing the case of the HM cluster. The outermost population is typically dispersed by the very first close stellar encounter in each simulation, giving rise to a continuous distribution of semi-major axes reaching beyond 1 pc. Such wide orbits are unstable in the current Galactic environment and even more so in a dense cluster. The comets are rapidly lost by the two energy-perturbing agents identified above, that is, the stellar encounters and the cluster tide.

After some 200 Myr, the primordial populations have dispersed sufficiently that the comet orbits appear rather smoothly distributed, albeit retaining a broad central peak covering a ~ 2000–15 000 au. The subsequent time evolution sees the outer population, a> 10 000 au, diminish at a rate similar to the central population. This means that a steady state is reached, whereby this outer population remains as a transit stage of comets migrating from the inner parts and eventually departing from the solar system due to energy perturbations. The semi-major axis distribution for a> 20 000 au is seen to evolve toward a power-law slope of − 2 corresponding to a flat energy distribution, which is typical of a diffusion process with an absorbing wall near 1 /a = 0 caused by our definition of departures.

The innermost part of the comet cloud, a< 5000 au, is rapidly populated during the first 50 Myr. The core of this population, a ≲ 3000 au, then remains essentially inert, being depopulated only in those individual simulations where very strong and close encounters occur. This core appears to form as a rule rather than an exception, comprising roughly 20% of the surviving comets from the primordial inner cloud, while the primordial intermediate and outer clouds each contribute one and two orders of magnitude fewer comets (see below). The core is thus more than twice as populated as the outward migrators beyond 20 000 au at the time (T = 400 Myr), when we stop the integrations.

From the plotted results, we also have an indication of what would happen if we had continued the integration further in time. The peak would remain close to 5000 au, and the curve would flatten out at smaller semi-major axes, while it would continue to be shifted downward at larger semi-major axes. Thus the predominance of the quasi-inert core would become further accentuated, as the total population of the cloud continues to decrease.

The evolution of the cloud in the IM cluster case – shown in the right panel – is basically similar, but some differences are easily seen. After 50 Myr the structure undergoes very small changes. As a rule, the Sun has then left the cluster or migrated out of the central region for most of the time. We noted above that this leads to a fossilized structure of the cloud, which we here see represented by the histograms in red, blue and green. The inner core is less pronounced than in the HM case. The outermost part of the cloud is cut at a ~ 50 000 au, since there are no more perturbations large enough to replenish these orbits from the inside.

Fig. 8 Mean semi-major axis distributions of the surviving comets at the end of the integrations, combined from all simulations, for the three primordial populations as in Fig. 6, using the same colours and styles of the curves. Vertical dashed lines represent the initial values. Differential distributions are shown by histograms, using a common normalisation to the size of the initial populations. The curves are shaped by energy perturbations caused by the external agents. Statistical noise is seen mainly for the inner and intermediate primordial populations. Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster.

We illustrate in Fig. 8 the separate mean semi-major axis distributions for the surviving comets of the three primordial populations. The HM cluster is shown to the left and the IM cluster to the right. The median values of the survivors in the HM cluster are less than the initial values – 4510, 8000 and 15 700 au for the inner, intermediate and outer populations, respectively. These shifts are due to the preferential loss of comets reaching large semi-major axes and the relative safety of comets diffusing toward smaller values. The shifts are smaller in the IM cluster, especially for the inner and central populations. This is likely due to the absence, in most cases, of a long-term energy diffusion.

Thanks to the common normalization used, the fact that the outer population has the smallest number of survivors in both clusters is clearly displayed in Fig. 8. The mean survival probabilities in the HM cluster are 2.1, 0.6, and 0.1% for the three respective populations. In the IM cluster these values are 11.7, 4.1, and 0.8%. All these values are significantly larger than the corresponding medians (Tables 1 and 2), because there is significant spread between the results of different simulations, and the survivors are concentrated to the minority that had the smallest external effects.

The sums of the three differential distributions yield the distributions shown at T = 400 Myr in Fig. 7. Each of these has a roughly triangular shape in the log-log diagrams, with a maximum at the initial value of the semi-major axis. For the HM cluster we see a steeper slope for larger than for smaller values. As noted above, the steeper slope is close to − 2, and the core population created by inward migration has contributions differing by roughly one order of magnitude between the inner, intermediate, and outer primordial populations. For the IM cluster the distributions are more symmetric around the maximum until the cut at large semi-major axes is reached. The slopes are higher on both sides of the maximum than in the HM cluster case.

Fig. 9 Statistics for the IM cluster over the distribution of simulations (grey-shaded histogram) and the survival percentages of comets. The comet populations are shown both combined (thick black curve) and for each population (coloured curves). The survival statistics are computed as the arithmetic mean using a boxcar, and then smoothed by a Gaussian kernel. Left panel: statistics as a function of the time when the Sun is ejected from the cluster, using a boxcar width of 50 Myr and Gaussian kernel σ = 10 Myr. Cases where the Sun is not ejected from the cluster are shown to the right of the vertical dashed line. Right panel: statistics as a function of the Sun’s nearest distance to the cluster centre, computed in logarithmic bins using a boxcar width of 0.3 dex and Gaussian kernel σ = 0.05 dex.

We have already noted a few special features of the evolution experienced within the IM cluster. One is that the Sun tends to leave the cluster during our simulations, whereby the influence of the cluster on the Oort Cloud is terminated. The other is that the IM cluster is equipped with a high-density central core, which acts as a very efficient pitfall to Oort Cloud comets, in case the solar orbit enters into its vicinity. To explore the influence of these features on the survival of the Oort Cloud, in Fig. 9 we illustrate the relevant statistical properties: histogram distributions of the time when the Sun is ejected from the cluster and the minimum periapsis distance of the solar orbit, and together with these the variations of the Oort Cloud average survival probability with the parameters in question.

The most striking fact revealed by the two panels of Fig. 9 is that the minimum periapsis distance effectively governs the Oort Cloud survival probability, while the time of ejection of the Sun does not exhibit any similar influence. It is clear that any approach of the Sun to less than 0.2−0.3 pc from the cluster centre leads to the loss of almost the entire Oort Cloud including the comets from the primordial inner cloud. On the other hand, for periapsis distances approaching 1 pc most of the inner and intermediate parts will survive. We interpret this to reflect the tidal action of the cluster core, whose radius is approximately 0.1−0.2 pc (see Fig. 2). Since the median of the minimum periapsis distance is seen to be close to 0.1 pc, such tidal losses can fully explain the generally low survival probability of Oort Cloud comets in the IM cluster.

Of course, when the Sun penetrates into or close to the cluster core, it can also experience close stellar encounters that strip comets away from the Oort Cloud. We showed in Table 2 that the median survival probability is very low, both if we turn off the cluster tide, and if we turn off the stellar encounters. This demonstrates that stellar encounters do have an influence. However, the tidal effect is probably the dominant effect. In our model, only one star will be closely encountered during each periapsis passage of the Sun, and it seems unlikely that the orbit with the smallest distance from the cluster centre will invariably involve an encounter that is efficient enough to strip away almost the whole cometary cloud. Thus, some scatter would be expected, making the survival curves less smooth and monotonic if stellar encounters were dominant.

The wavy pattern exhibited by the survival probability with respect to time of ejection reflects the limited statistics. Although the curves are smoothed, their maxima are strongly influenced by the occasional simulations, where the minimum periapsis distances are large and many comets survive. Ejection times close to that of one of these simulations will generate a relatively large survival probability in the calculated distribution. The waves for different initial cloud populations are in phase, because these groups are simulated with identical solar orbits and encounters. The absence of any systematic decrease of the survivability with time of ejection merely shows that the Sun may spend a long time in the cluster with its Oort Cloud intact, provided that it does not penetrate into the cluster core.

3.4. The solar orbit

Fig. 10 Distributions of initial and final solar orbits for all our simulations using the full dynamical model. Pericentre (q) and apocentre (Q) distances of the rosette orbits are plotted on the axes of each diagram. Left panels: the high-mass cluster. Right panels: the intermediate-mass cluster. Top panels: initial orbits. For the HM cluster, these orbits start from the half-mass radius (r = 4.87 pc) with some random radial velocity, so we always have q< 4.87 pc <Q. The red crosses denote orbits, where the Sun was ejected from the cluster during the simulation, while the rest are denoted by black diamonds. Bottom panels: final orbits. Only orbits where the Sun survived as a cluster member are shown.

Let us finally present some results on the evolution of the solar orbit in the cluster, caused by the impulses received from the encountering stars. Figure 10 shows two scatter diagrams of apocentre versus pericentre distances at the beginning of our integrations (upper panels) and at the end (lower panels). The left pair of panels refers to the HM cluster and the right pair illustrates the IM cluster. Each symbol represents one of the 1000 simulations with the full dynamical model. The red symbols in the initial distributions mark those solar orbits that were not stable, leading to the Sun being ejected from the cluster before the end of the integrations (this fate was registered, if the Sun moved beyond the tidal radius, given in Table A.1, in each cluster).

The cases of ejection amount to 5.8% of the simulations for the long-lived, HM cluster. This verifies the expectation that the Sun, as a relatively massive cluster star, suffers only a small risk of ejection within 400 Myr. For the IM cluster the situation is the opposite. The cluster as a whole dissolves on a much shorter timescale and the Sun is no exceptional star. We find that the median survival time of the Sun as a cluster member is about 100 Myr, while in 80% of the simulations the Sun is ejected within 200 Myr. After 400 Myr, the Sun remains in only 4.8% of the simulations.

For both cluster models, the initial and final distributions of solar orbits differ markedly from each other. Most of the time the Sun is pushed outwards in the cluster, since both apocentre and pericentre distances show increasing trends – more pronounced for the apocentre distance. For a minority of cases in the HM cluster, the solar orbit evolves into a smaller apocentre distance or a lower effective eccentricity. In the IM cluster, most of the remaining solar orbits stay close to the cluster centre.

We caution that the trend toward the outskirts of the HM cluster may be affected by our manner of selecting the stellar encounters. Choosing each time the encounter with the largest strength parameter (Sect. 2.3) means that the massive stars are favoured as encounter partners. Statistically, during binary encounters, energy per unit mass flows from the more massive to the less massive partner – a basic reason for mass segregation causing the concentration of high-mass stars to the cluster core and the preferential escape of low-mass stars (Spitzer 1969). Had we included all the encounters, in which case the simulations would have been more realistic, the outward migrating trend for the Sun would likely have been balanced by many interactions with low-mass stars and thus reduced.

Another word of caution is justified concerning the IM cluster. Here, again, our result may be biased by our modelling of the cluster. In reality, the central region of the cluster becomes enriched in massive stars, but we neglected this mass segregation. Thus, our encounter selection – while preferring massive partners – did not favour these as much as it should have done in a realistic modelling. This may have made it easier for the Sun to stay in the central region rather than be expelled from it.

4. Discussion

4.1. Modelling issues

We have chosen to base one of our cluster models on M 67, which according to Hurley et al. (2005) started out with more than 20 000 stars. This was done in spite of the conclusion by Adams (2010) that a birth cluster with more than about 10 000 stars would threaten the stability of the planetary orbits. One reason not to worry is that Adams considered giant planets on their current orbits, while the Nice Model holds that the orbits of the giant planets had much smaller semi-major axes during the early epochs of solar system history – both in its original form (Tsiganis et al. 2005) and in later versions. This would obviously reduce their vulnerability to external perturbations. However, the intricate resonant clockwork of the Nice Model (Levison et al. 2011) might be upset by stellar encounters far more distant than those previously considered as disastrous. No analysis of this problem has yet been made to our knowledge.

We did not aim to survey the full spectrum of cluster sizes but took two examples of relatively rich systems that are statistically likely birthplaces of solar type stars (Lada & Lada 2003). We found that the intermediate-mass cluster is in a certain sense as hostile to the survival of a primordial Oort Cloud as the high-mass one. However, it is dangerous to extrapolate this to even smaller birth clusters. At some point, the cluster dissolves so quickly that the destructive influence on the Oort Cloud disappears. At which initial cluster mass this transition occurs remains to be found out.

Meanwhile, it is worth noting that a surviving Oort Cloud formed in a very dense cluster environment should be very tight, including semi-major axes far smaller than those that we have investigated (Fernández & Brunini 2000; Brasser et al. 2006). We have seen that such an inner core does survive in the clusters we modelled. An important issue is then how to activate this core and repopulate the outer halo after the star leaves the cluster. In a sense, a fossilized inner core would be ineffective: it would not provide any observable comets. We shall return to this discussion in Sect. 4.4.

As for the Sun, an active Oort Cloud was formed in the simulations by Levison et al. (2010), who considered a very small and short-lived birth cluster. In fact, judging from Adams (2010), many of the relevant cluster sizes would not be in absolute conflict with the nucleosynthesis requirement, since they would yield a random probability of ~10% for a relevant supernova explosion.

Another issue concerns the realism of our dynamical model for the external perturbations suffered by the Oort Cloud comets. Of course, only a full N-body simulation of the cluster may be considered fully realistic, but this is beyond the scope of our preliminary study. The synthetic model that we developed represents the cluster by its smooth potential field plus a sequence of two-body encounters within this field, which the Sun experiences with individual cluster stars. This is quick and efficient, but it certainly departs from reality. We have assumed the Sun and the comets to experience, in addition to the smooth cluster potential, the fluctuating component caused by just one passing star at a time. In reality there are typically several nearby stars contributing to this fluctuating field at the same time. Hence, concerning the loss of comets, in many cases the fate of a comet may critically depend on the details of the tidal field, so that our approximation may either save comets from leaving the Sun or stimulate their escape, depending on the circumstances.

Although it does not seem likely that any serious systematic errors occur as a result of our simplifications, a full answer cannot be found except through time-consuming N-body simulations. One technique would be to relax the blocking of overlapping encounters to see, if by allowing up to 5–10 simultaneous encounters and treating these N-body systems accurately, we would get statistically different results on the loss of comets. However, this would mean an additional major effort, which still would not solve all problems.

4.2. Simplifying assumptions

The most important simplification that we used may be the assumption of a step function for the initial population of the Oort Cloud. In fact, the Oort Cloud was not built instantly. It is clear from the works of Kaib & Quinn (2008) and Brasser & Morbidelli (2013) that the emplacement of comets from the planetary region into the Oort Cloud took several hundred Myr. Thus, what we call a primordial Oort Cloud should have been enriched in new members for a time comparable to the cluster lifetime or the full length of our simulations.

Connected to this is another assumption, namely, that the Sun still resides in its birth cluster at the starting time that we use, 100 Myr after the formation of the Sun and the cluster. For the HM cluster, the risk of ejection of the Sun during this early interval is negligible, but not so for the IM cluster. This has an initial mass of 820 ℳ_⊙, and at the starting time the mass has decreased to 710 ℳ_⊙. This decrease amounts to 13%, which we take as a rough estimate of the risk of early ejection of the Sun.

Hence, in the IM cluster case, there is a 13% chance that some of the comets transferred into the Oort Cloud during the first 100 Myr as well as all those transferred at later times would not feel any effects of the birth cluster. However, the comets in question – likely being the majority of all Oort Cloud comets – would then be emplaced without the help of the birth cluster. The Oort Cloud would then likely have much less of an inner core than if the Sun had stayed in the cluster, and the creation of Sedna-type objects might be strongly curtailed.

With the complementary probability of 87%, we have underestimated the cluster influence on the early emplaced comets and yet strongly overestimated the influence on those that were emplaced at later times. It seems clear that the overestimates dominate as an error source. Again, one would have to distinguish between two parts of the Oort Cloud – the comets that were emplaced inside the birth cluster, which did experience its destructive effects, and those that were emplaced after the Sun had left, which did not benefit from the cluster in populating the inner core.

At any rate, we conclude that the Sun could not spend a longer time than approximately 50 Myr in an IM cluster with an already formed Oort Cloud left intact except under special circumstances. Such circumstances would include a solar orbit that kept the Sun and the Oort Cloud constantly outside the central part of the cluster.

Finally, we stress that there are three categories of birth clusters, which we have estimated to be approximately equally likely for the Sun in terms of the initial cluster mass function. The LM case is, however, less probable in view of the nucleosynthetic evidence for an early supernova in the neighbourhood (see the discussion above). In the HM case our model should be reasonably good; the IM case is discussed immediately above; and in the LM case the destructive effects of the birth cluster would likely be much smaller.

In fact, we have made a few additional approximations that likely caused us to underestimate the losses of comets, thus yielding conservative estimates. The first is that we neglected mass segregation in our model cluster. Hence we downplayed the risk for the Sun encountering a massive star in the high-density environment near pericentre, which would have had dire consequences for the entire Oort Cloud.

The second is that we neglected binary stars. M 67 is known to be rich in binaries (Richer et al. 1998), and thus other similar clusters may be suspected to be similarly binary-rich. Binary stellar systems are also known from dynamical simulations to form and dissolve within star clusters including clusters of much smaller masses (Giersz & Heggie 1997). By ignoring binaries, we have artificially increased the number of potential encounter partners of the low-mass type, while entirely neglecting a kind of partner that would have had a great capability to transfer energy and momentum to the Sun’s motion – and, similarly, to destroy the Oort Cloud. In particular, “hard” binaries in close orbits are known to statistically give energy away to the encounter partner, thus providing an energy source for the dynamical evolution of the cluster, as reviewed by Elson et al. (1987).

A third approximation is that we treat only a rather small number of subsequent encounters in each simulation. This number might be increased without introducing overlapping encounters, but a more fundamental issue is that we select the strongest encounters in terms of the impulse imparted to the Sun per unit mass within the approximation of the classical impulse approximation. Certainly, all kinds of encounters may occur, but we systematically disfavour the weaker ones involving the less massive partners. As already remarked, this creates an exaggerated trend for the Sun to move outward in the cluster. Hence, statistically, the Sun spends too much time in the outer, less populated regions, so we underestimate the risk of the strong perturbations that characterise the inner parts. Thus, our treatment should underestimate the total number of lost comets from the cloud in that an overly small number of stellar interactions are considered, which in fact are not necessarily the most effectively destructive ones. Our estimated destruction rates for the Oort Cloud are thus conservative in this respect.

Yet another neglected phenomenon could have increased the number of remaining comets. Some of the lost comets may be picked up by the encountering stars, and this number would likely increase, if we would treat overlapping encounters. Even though most lost comets rather become cluster vagabonds, in due time these could also be picked up by some cluster star, including the Sun, in a way similar to the formation of stellar binaries in clusters.

4.3. Extra-solar Oort Clouds

All this begs the question of whether or not by assuming that the Sun had a primordial Oort Cloud, one should assume that other cluster members of similar types were also equipped with such primordial cometary clouds. If so, the sort of cluster environment that we consider here might stimulate a certain exchange of comets between different stars (Zheng et al. 1990; Levison et al. 2010). We cannot say how efficient this process would be using our simplified model. On the other hand, it is clear that picked-up Oort Cloud members would be particularly vulnerable to being lost during following encounters, since these would occupy relatively loosely bound orbits.

Concerning the question, if Oort Clouds may be a characteristic feature of Galactic disk stars, we note that Stern et al. (1991) searched for extra-solar Oort Clouds around 17 nearby stars by looking for IR excess radiation using IRAS low-resolution data and an S/N-enhancing method, with negative results which, however, may be ascribed to the limited sensitivity. Black (2010) has extended this search on the basis of the IRAS sky survey to all F and G dwarfs, augmented with all other stars (then) known to have planetary systems, within 50 pc distance from the Sun. While challenging, since the sensitivity of the IRAS is again a severely limiting factor, no positive identifications of Oort Clouds around other stars were reported.

Some cold, dusty outer disks have been found and studied by means of the Herschel telescope, around young stars, also with planetary systems. One example is the A5 V star HR 8799 with four known planets, at a distance of 40 pc and an age estimated at 20 to 50 Myr, which has a central warm dust component, an outer cold component extending from 90 to 300 au and an outer component of small grains extending beyond 1000 au (Matthews et al. 2014). The evidence for any clumping in this halo is, however, meagre. Most observations of debris disks around young stars are limited to A-type stars and to rather small radial distances from the star. One interesting example is, however, the F5/6-type star HD 181327, a member of the 20 Myr old β Pictoris moving group, for which ALMA observations disclose a ring-like CO gas disk, in addition to the dust ring, with a halo extending out to 200 au, and a CO + CO₂ cometary composition (Marino et al. 2016).

Neutron stars capturing comets when passing through the Oort Clouds of other stars have been suggested to provide indications on the (non-)existence of extra-solar cometary clouds. Shull & Stern (1995) thus proposed that such events should generate repeating bursts of soft gamma-rays, and estimated that the absence of such Galactic events indicated that, at most, a few percent of the Galactic stars have Oort Clouds. It is, however, very questionable whether accretion of comets onto neutron stars would occur abruptly enough to generate such bursts; weaker emission, more extended in time and at lower energies, seems more probable (we thank Dr. J. Poutanen for making this point). We conclude that presently no observational limits may as yet be set on how frequent Oort Clouds are around stars.

4.4. Concluding remarks

Our basic result is that, if the Sun was born as a member of a relatively rich stellar cluster, a significant part of a primordial Oort Cloud would not likely survive the time in the cluster until the Late Heavy Bombardment. The relative extent of this depletion depends on the detailed orbital evolution of the Sun within the cluster. In the case of a massive depletion, the formation of the present Oort Cloud as a consequence of a late planetary migration within the Nice Model (Levison et al. 2011) appears to be a viable scenario, provided that the Sun had then already left the cluster, or was about to do so. Such a scenario has been explored by Brasser & Morbidelli (2013).

In a very massive cluster the escape of primordial Oort Cloud comets is mainly caused by the disrupting effect of stellar encounters. The eccentricity pumping due to the cluster tide plays a role only for a minority of comets in the inner part of the cloud. In the outer part the energy perturbations caused by the tide may constitute an important source of comet losses. In clusters of lower mass the latter type of tidal perturbations provide a major loss mechanism in all parts of the Oort Cloud. If the cluster mass is very high, the number of comets penetrating to within 5 au of the Sun amounts to ≃4% of the primordial inner cloud members. For the lower cluster mass, this fraction decreases to about 2%. If the primordial Oort Cloud had a high mass, this might cause an important cometary bombardment of the terrestrial planets (including a late veneer of H₂O) well before the time of the LHB.

The influence of the stellar encounters on the solar orbit can be seen in Fig. 10. From this and the solar escape rate, we conclude that the Sun typically receives a cumulative impulse of several km s^-1. Since the number of encounters per simulation is only about 20, in a random walk there must be some individual kicks experienced by the Sun that amount to about 1 km s^-1, thus strong enough to unlink most of the Oort Cloud comets. The solar impact parameters of those stars are typically much less than the Sun-comet distances. Thus, the main mechanism of comet escape in our model is that the Sun is kicked away from its comet cloud rather than individual comets being kicked away from the Sun.

By restricting our simulations to a small number of stellar encounters, we may have introduced too large a statistical dispersion of the results. In particular, the simulations that left the Oort Cloud – particularly for the intermediate-mass cluster – with much more comets than the median can be regarded as chance selections of solar orbits with weaker tidal effects and encounter sequences than normal. We have shown that the minimum periapsis distance q_min of the Sun plays a decisive role for the survival frequency in the Oort Cloud – the larger q_min, the more survivors. According to our results, the category of outcomes with the larger survival rate makes a significant though not dominant contribution, but some caution is warranted, especially concerning the results for the IM case, until more realistic cluster simulations can be made. Tentatively, there is no reason to suspect that we would have exaggerated the comet-loss rate in the IM cluster by selecting initial solar orbits with overly small q_min. The selection was made on the basis of the general density profile of the cluster, so we neglected the fact that the Sun, being a relatively high-mass star, would tend to prefer the central region at the time we started our simulations and the cluster was already considerably relaxed. Again, however, only fully realistic simulations including binarity and other special phenomena will provide the full answer to this question.

A primordial Oort Cloud may consist of comets that originated in the accretion zone of the giant planets. However, if the planetary orbital instability of the Nice Model actually happened very early, as was recently argued by Kaib & Chambers (2016), the Oort Cloud originating from the trans-planetary disk would also be primordial. It may have been formed after the birth cluster was dispersed, if this cluster contained only about a hundred stars. Such a birth cluster may not be totally excluded by the argument of nucleosynthesis providing the short-lived radionuclides for the solar nebula, even though somewhat special circumstances would be called for in the case of the Sun. If so, the formation scenario of Brasser & Morbidelli (2013) would be a relevant model except that the timing of the event is different.

On the other hand, if the birth cluster was of the IM or HM kind, we have shown that the primordial Oort Cloud would largely survive only as a tight, inner core. Models of Oort Cloud formation in such a dense stellar environment indeed predict an initial cloud structure dominated by such a core (Brasser et al. 2006). Therefore, in such a case the existence of the present outer halo, from where the observed comets can be transferred by the Galactic tide, requires a mechanism of energy transfer that can activate the core from its inert state. The alternative would be a late planetary instability as investigated by Brasser & Morbidelli (2013), in which case it is reasonable to assume that the Sun had already left the birth cluster.

To be specific, taking the Galactic disk tide to be the mechanism for bringing Oort Cloud comets into the inner solar system, the tidal torquing timescale (Heisler & Tremaine 1986) is found to be longer than the age of the solar system, unless a ≳ 5000 au. Stellar perturbations would not change this result drastically. We find that such comets in general would not survive the dwelling time within the birth cluster. A fossilised, primordial inner core of comets with a ≲ 3000 au might exist, possibly including the Sedna-type objects, but it would not be able to produce the observed, new Oort Cloud comets. The gap in orbital energy would have to be bridged by an as yet unidentified mechanism.

Returning to the issue of extra-solar Oort Clouds, our results suggest that such clouds would often be severely truncated by the effects of the birth clusters and therefore be much smaller than usually thought. Quite likely, Oort Clouds in general would mostly stem from trans-planetary planetesimal disks, and thus, their existence depends on the way extra-solar planets have migrated.

Acknowledgments

This work was supported by the Polish National Science Center under Grant No. 2011/01/B/ST9/05442 as well as Grant No. 74/10:2 of the Swedish National Space Board. T.N. acknowledges support from Swedish National Space Board (Rymdstyrelsen) and Anna & Allan Löfbergs stiftelse. Our thanks go to Ross Church, Mirek Giersz, Alessandro Morbidelli and Juri Poutanen for enlightening discussions, and to the referee, Nathan Kaib, for raising points of crucial importance for the interpretation of our results. Computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX), Lund University Center for scientific and technical computing (LUNARC), and High Performance Computing Center North (HPC2N), under projects s02011-10, snic2013-11-24, and SNIC2015-1-309.

References

Appendix A: Calculation of stellar cluster structure

We prescribe that the model cluster should have a given mass ℳ_cl and be situated in the Galactic potential at distance r_G from the Galactic centre. This allows the calculation of a tidal radius of the cluster, R_t. Using a limiting energy E_t, taken to be the energy of a circular orbit situated at a distance from the cluster centre r = R_t, our King model is calculated from a distribution function ϕ, which can be expressed as $\begin{array}{ccc}\mathit{\varphi }\mathrm{\left(}\mathit{E}\mathrm{\right)}& \mathrm{=}& \end{array}$(A.1)where E = U + v²/ 2 is the total energy of a test particle, and a and b are positive constants of the model. We rewrite this distribution function when E<E_t as a function of velocity and radial position, $\begin{array}{ccc}\mathit{\varphi }\mathrm{\left(}\mathit{v,r}\mathrm{\right)}& \mathrm{=}& \end{array}$(A.2)and solve for the density ρ, $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{v}}_{\max }\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{\varphi }\mathrm{\left(}\mathit{v,r}\mathrm{\right)}{\mathit{v}}^{\mathrm{2}}\mathrm{d}\mathit{v}\\ & \mathrm{=}& \mathrm{4}\mathit{\pi a}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{v}}_{\max }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{v}}^{\mathrm{2}}\left\{{\mathrm{e}}^{\frac{\mathit{b}}{\mathrm{2}}{\mathit{v}}_{\max }^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{b}}{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right\}\mathrm{d}\mathit{v}\\ & \mathrm{=}& \end{array}$(A.3)

The velocity distribution can be evaluated by multiplying (A.2)by 4πv², giving locally a truncated Maxwellian distribution. The structure of the cluster, its distribution of mass and the gravitational potential, is finally computed by solving the integral (A.3)simultaneously with Poisson’s equation, $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{U}}{\mathrm{d}{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r}}\frac{\mathrm{d}\mathit{U}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{4}\mathit{\pi G\rho ,}& & \end{array}$(A.4)using the parameter substitutions suggested by King (1966, Eqs. (19)–(22)). The equations are solved from the cluster centre, where the potential energy U₀ is taken as a free parameter. The free parameters a,b,U₀ thus represent a reformulation of the classical free parameters W₀, r₀ and ℳ_cl required to fit a King model to an idealized stellar cluster.

The input parameters and resulting properties of the optimal model are given in Table A.1.

Appendix B: Implementation of stellar encounters

Each list covers a time interval of 40 Myr, and it comes from a preliminary integration of the solar orbit in steps of 3 kyr, yielding the mean encounter frequency, ⟨ n ⟨ v_rel ⟩ ⟩, where n and ⟨ v_rel ⟩ are interpolated from the solution shown in Fig. 4. For each step, we let an encounter occur within b_max = 20 000 au, if $\begin{array}{ccc}{\mathit{b}}_{\max }^{\mathrm{2}}\mathit{\pi }{\mathrm{\int }}_{{\mathit{t}}_{\mathit{i}}}^{{\mathit{t}}_{\mathit{i}\mathrm{+}\mathrm{1}}}\mathit{n}\mathrm{⟨}{\mathit{v}}_{\mathrm{rel}}\mathrm{⟩}\mathrm{d}\mathit{t}\mathit{>}\mathit{\xi ,}& & \end{array}$(B.1)where t_{i + 1} = t_i + 3 kyr, and ξ ∈ [0,1] is a random number drawn for each step from a uniform distribution. This typically generates about 100 encounters. For each such encounter, we generate a random stellar mass from the IMF as described in Sect. 2.1.3, evolving it to the age of the cluster at the current time. If the star is found to have evolved through a supernova phase, it is discarded and a new value is drawn from the IMF and again evolved to the current time. The impact parameter for each selected encounter is determined randomly using $\mathit{b}\mathrm{=}{\mathit{b}}_{\max }\sqrt{\mathit{\xi }}$, for a new random ξ ∈ [0,1]. Finally, the star is assigned a velocity v_rel relative to the Sun as in Sect. 2.2.

After this preliminary selection, we know the times and parameters of all the potential, upcoming stellar encounters, and we are able to pick one of them based on the strength parameter S = ℳ/(v_relb). When modelling the encounter in a two-body scattering problem, the approximation would be to let the encountering star aim from infinity at a position on a heliocentric circle with radius b in a plane perpendicular to the direction of approach (called the impact plane). However, under the influence of the cluster potential, this straight-line approximation cannot be used. Instead, we realise the closest approach by first choosing randomly the direction of relative motion, which defines the impact plane, and then placing the star on the impact plane at a distance b from the Sun at a random azimuth.

The orbits of both stars are then integrated backward in time in the cluster potential with no mutual gravitational interaction, until the distance between the stars is d_start. The time and geometric configuration at this stage are stored as the initial state of that encounter. In the few cases of very slow encounters, where the backward integration overlaps the previous encounter, the setup is considered as failed and is discarded. We then select the second largest value of S and repeat the calculation of the initial state.

All Tables

All Figures

	Fig. 1 The number of stars (solid, black lines) in our high-mass and intermediate-mass cluster models, and their half-mass radii (dashed, red lines), versus time from cluster formation. A logarithmic scale is used for the number of stars and a linear scale for the radius. The time has been rescaled for the IM cluster to take the expected influence of GMCs into consideration (see text).
In the text

Fig. 2 Density distribution in the template clusters. The dashed black curve indicates the distribution of the high-mass cluster, while the solid coloured curves represent time steps in the intermediate-mass model, 100 Myr apart, with the average time of each step shown in the legend. The central mass density (left end of the plot) decreases monotonically with age of the model. The horizontal dash-dotted line indicates the Galactic mid-plane density for comparison.

In the text

Fig. 3 Distribution of relative velocities in our encounter model for the stellar cluster. The five curves in the left panel represent, from top to bottom, the maximum allowed relative velocity (solid, black), its 90th percentile (dotted, red), mean (solid, blue) and 10th percentile values (dotted, black). Shown as vertical (black) lines are the core (dotted) and half-mass (dashed) radii. Left panel: relative velocity distributions for the high-mass cluster. Right panel: relative velocity distributions for the intermediate-mass cluster at two time steps. The first step represents the time interval 100–125 Myr. The second time step represents the time interval 300–325 Myr, and is shown only by a thick dashed black line representing the maximum relative velocity, a thick dashed blue line representing the mean, and a vertical thick dashed line representing the half-mass radius at this time. The core radius is not shown at this time, but is available for all time steps in Table A.1.

In the text

Fig. 4 Encounter frequency versus distance from the cluster centre, shown for the high-mass cluster (dashed line) and selected time steps for the low-mass cluster (coloured solid lines). The encounter flux at the centre of the cluster decreases with increasing age for the intermediate-mass cluster – colours are the same as in Fig. 2. The horizontal dash-dotted line illustrates for comparison the current encounter frequency of the Sun with Galactic field stars (Rickman et al. 2008).

In the text

Fig. 5 Solar orbit during the first 80 Myr in a simulation of the high-mass cluster, with grey circles indicating the core and half-mass radii. The solar orbit is illustrated by a line that switches colour (from violet to blue, green and yellow) every time a star is encountered. The starting points of the interloping stars are numbered and indicated by a star symbol, while end points are indicated by x. The starting point of the solar orbit is indicated by a circle of radius 20 000 au. The solar orbit is by construction initially confined to the x-z plane, and departs significantly from that plane only after a close interaction with star number 3.

In the text

	Fig. 6 Distribution of comet survival probability for the three different bins of initial semimajor axis, ao = 5000 au (red), 10 000 au (blue), and 20 000 au (orange). Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster.
In the text

Fig. 7 Time evolution of the combined comet population of all simulations with the full dynamical model. The fractional populations are counted with respect to the total initial population. Results were sampled at 10 Myr timesteps, using 100 uniform logarithmic bins for a ∈ [500,200 000] au. The locations of the three initial sub-populations are indicated by vertical dashed lines. Here, snapshots are shown at times T = 10, 50, 200, and 400 Myr, from top to bottom, using black, red, blue and green colours, respectively, where the final time step is also shaded in grey. Comets beyond the plot limits have been summed at the edges. Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster. The three remaining peaks in the final distribution reflect the fact that the vast majority of the stars in this cluster model have by then left the cluster with fossilised structures of their cometary clouds.

In the text

Fig. 8 Mean semi-major axis distributions of the surviving comets at the end of the integrations, combined from all simulations, for the three primordial populations as in Fig. 6, using the same colours and styles of the curves. Vertical dashed lines represent the initial values. Differential distributions are shown by histograms, using a common normalisation to the size of the initial populations. The curves are shaped by energy perturbations caused by the external agents. Statistical noise is seen mainly for the inner and intermediate primordial populations. Left panel: the high-mass cluster. Right panel: the intermediate-mass cluster.

In the text

Fig. 9 Statistics for the IM cluster over the distribution of simulations (grey-shaded histogram) and the survival percentages of comets. The comet populations are shown both combined (thick black curve) and for each population (coloured curves). The survival statistics are computed as the arithmetic mean using a boxcar, and then smoothed by a Gaussian kernel. Left panel: statistics as a function of the time when the Sun is ejected from the cluster, using a boxcar width of 50 Myr and Gaussian kernel σ = 10 Myr. Cases where the Sun is not ejected from the cluster are shown to the right of the vertical dashed line. Right panel: statistics as a function of the Sun’s nearest distance to the cluster centre, computed in logarithmic bins using a boxcar width of 0.3 dex and Gaussian kernel σ = 0.05 dex.

In the text

Fig. 10 Distributions of initial and final solar orbits for all our simulations using the full dynamical model. Pericentre (q) and apocentre (Q) distances of the rosette orbits are plotted on the axes of each diagram. Left panels: the high-mass cluster. Right panels: the intermediate-mass cluster. Top panels: initial orbits. For the HM cluster, these orbits start from the half-mass radius (r = 4.87 pc) with some random radial velocity, so we always have q< 4.87 pc <Q. The red crosses denote orbits, where the Sun was ejected from the cluster during the simulation, while the rest are denoted by black diamonds. Bottom panels: final orbits. Only orbits where the Sun survived as a cluster member are shown.

In the text