© ESO 2021

1. Introduction

The Sun is orbited by a cloud of minor bodies in a spherical distribution between ∼10⁴ au and the Hill radius of the Solar System in the Galactic potential. This cloud of objects, generally referred to as the Oort (1950) cloud, may not be unique to the Solar System: Other stars may also have their own ‘Oort’ clouds (Stern 1989; Torres et al. 2020; Torres 2020). The same circumstances that lead to the Sun’s Oort cloud could exist for many nearby stars, in which case Oort clouds should be rather common.

Crucial to forming an Oort-like cloud is the presence of a debris disk with relatively massive objects (which we call planets; Brasser & Morbidelli 2013). These planets perturb the low-mass minor bodies (which we call asteroids) and drive their eccentricities to high values (≫0.9) while preserving the closest approach distance to the parent star. Upon each orbit, the asteroid is kicked into a wider and more eccentric orbit while preserving pericenter distance. Eventually, the Galaxy’s tidal field, in which the planetary system orbits, starts to circularize the orbit of the minor body (Duncan et al. 1987). If the tidal field effectively reduces the orbital eccentricity of the minor body, it can remain bound to the parent star and become a member of the local Oort cloud. If the body becomes unbound, it continues as a free-floating asteroid, or sōlus lapis, in the Galactic potential (see Torres et al. 2019a, 2020).

The recent arrival of two such sōlī lapidēs, 1I/(2017 Q3) ’Oumuamua (Bacci et al. 2017; Meech et al. 2017a,b) and 2I/(2019 Q4) Borisov (2019) (and maybe even several others; Froncisz et al. 2020) in the Solar System provides an incentive to study the extent to which an evaporating Oort cloud or circumstellar debris disk of nearby stars contributes to the frequency and characteristics of such objects.

It has been suggested that the two observed sōlī lapidēs originate from the Oort clouds of other stars (Portegies Zwart et al. 2018a; Raymond et al. 2018; Moro-Martín 2019), were ejected from a proto-planetary disk (Moro-Martín 2018), from a binary system (Jackson et al. 2018), or from a stellar association (Feng & Jones 2018), or were launched as a tidal fragment (Zhang & Lin 2020). At least ’Oumuamua and Borisov do not seem to originate from the Sun’s Oort cloud (Meech et al. 2018; Higuchi & Kokubo 2020).

Several studies have tried to identify the host (Zhang 2018; Bailer-Jones et al. 2020), but despite the precision of the Gaia-DR2 astrometry (Gaia Collaboration 2018) and the accurate ephemera for the two sōlī lapidēs, efforts remain unfruitful. Another complication in tracing back the interloper to its host star is our lack of knowledge regarding the age of the object. Stars as well as sōlī lapidēs move around in the Galaxy and tracing back one to another can be arduous; much in the same way as it might be hard to identify the owner of a particular dog in a dog-walking area.

Nevertheless, there is hope in identifying the host of sōlus lapis when ejected asteroids linger around the source along its Galactic orbit (see also Torres et al. 2020 and Chap. 4 of Torres 2020). On the other hand, sōlī lapidēs remaining in the orbit of their host star may not appear to be launched from the star itself, but from some distance away. The asteroid may therefore have come from a somewhat different direction and with a different velocity from that of the host star, making it hard to connect the two. The orbital velocity of an asteroid in orbit around its parent star at the Hill radius (dv/v ≡ δv ≲ 0.1 km s⁻¹) is considerably smaller than the orbit of the star around the Galactic center (∼240 km s⁻¹). While the host star orbits in the Galactic potential, the low-velocity ejected asteroids will form extended leading and trailing arms along the orbit of their host. This is referred to as the proximity argument.

Rather than trying to identify the host star of a sōlus lapis, we use the proximity argument to determine the probable association between the host star and its sōlī lapidēs. We do this by simulating the distribution of asteroids that have become unbound from their parent star. Although these asteroidal tails may be quite extended (several kiloparsecs), we limit ourselves to the 200 stars nearest to the Sun (within about 16 pc). Torres et al. (2020) performed a similar study for the nearest half dozen stars over a time frame of 20 Myr. Here we extend this latter study to the nearest 200 stars. We assume that these stars had a disk of planets and debris sometime in the past. In the simulations, Oort clouds form through planets that eject asteroids. Most of these become unbound and orbit the Galaxy. We subsequently study the Galactic distribution of these ejected asteroids and investigate the probability of them interacting closely with the Solar System.

We study the phase-space distribution of minor bodies from the moment they enter the conveyor belt and are kicked out beyond the vicinity of the local host star while sojourning the Milky way Galaxy. Our calculations start by integrating the 200 nearby stars backward in time in the Galactic potential for 1 Gyr. Once calculated backward, they are provided with a planetary system and a population of minor bodies as test particles in a disk. All stars with their planets and minor bodies are subsequently integrated forwards in time in the Galactic potential until the current epoch, where they are observed today, near the Sun. We call this the yo-yo scheme. During this orbital calculation we simulate the dynamical evolution of the planets and minor bodies. Once at the current epoch, the Oort cloud of each star has fully developed, and the majority of the minor bodies (and some of the planets) have been ejected on Galactic orbits. In this paper, we focus on the distribution of these sōlī lapidēs in the vicinity of the Sun. For safety, we continue the integration for another 100 Myr after the present time, in case the closest approach with the Sun happens to be in the near future.

2. Methods, initial conditions, and glossary

This numerical study focuses on a hypothetical population of interstellar asteroids that were liberated from their parent stars. After liberation, they float freely in the Galactic potential. Some of these sōlī lapidēs migrate to the Solar System or pass through the Oort cloud, but the majority will just hover around orbiting the Galactic center without ever interacting with the Solar System. The calculations depend on several assumptions, and these are reflected, in part, by the initial conditions chosen for this study. We realize that the reader may not be satisfied with some of the choices made, and for this purpose, all the run-scripts, data, and picture-generating scripts to reproduce this work are available through figshare and github. The source code for the AMUSE framework is available online¹ (see Portegies Zwart et al. 2018b).

2.1. Glossary of terms

The adopted terminology in this paper may be confusing at times. Therefore, we present a short glossary of terms used.

Asteroid: low-mass object (compared to the planets) in a bound orbit around a star.

Conveyor belt: Area in orbital parameter space (in particular semi-major axis and eccentricity) where an asteroid crosses the orbit of one or more of the major planets. These repeated perturbations induce small impulsive velocity kicks causing the orbits to drift to higher eccentricity and larger semi-major axis while preserving pericenter distance (sometimes called eccentricity pumping). The conveyor belt is indicated in Fig. 1 as the area between the two black curves.

Fig. 1. Orbital parameters of planets (blue studs) and 141864 asteroids (red points) around the 200 nearest stars 1 Gyr ago. Here the right-most light-blue curve gives the approximate inner edge of the Oort cloud using a mean local stellar mass of ∼1.0 M⊙ and a mean distance according to the Sun’s Hill radius in the Galatic potential (∼0.65 pc, red vertical line). The vertical dashed and dotted curves give the smallest and largest Hill radii experienced by any star throughout the simulations. The inner light-blue curve is given for a relative perturbation of δv = 10−8, and the outer light-blue curve gives δv = 10−5. The green curve gives the distance from the host star for which the circularization timescale of the Galactic tidal field is comparable to the eccentricity pumping timescale of a test particle in the conveyor belt. The orbital configuration at the present epoch, 1 Gyr after the start of the simulation, is presented here. This figure is generated with the Python script plot_semimajor_axis_vs_eccentricity_of_planets_and_comets.py. This script, the required data, and the scripts to generate any of the figures in this manuscript including all required data are available through figshare.

Frozen zone: area in parameter space (semi-major axis and eccentricity) where the orbits of asteroids are insignificantly perturbed by massive bodies (planets or stars) from shorter orbits, or stellar fly-bys.

Hill radius: We define the Hill radius as the distance to the star for which the gravitational force is dominated by the star. The Hill distance depends on the stellar mass and its orbit around the Galactic center. The Hill radius depends on direction, which results in an ellipsoid around the star. The force exerted by the star on objects within this volume exceeds the force from the rest of the Galaxy. The Roche radius (Kopal 1959, p. 136) is defined as the radius of a sphere with the same volume as the Hill ellipsoid.

Figure 1 presents the Hill radius as a solid red curve, here at roughly 2 × 10⁵ au. In our calculations, the Hill ellipsoid is only considered post-processing of the data. During the calculations, whether or not an asteroid (or planet) is bound to a star is calculated by integrating the equations of motion of the star, planets, and asteroids in the potential of the Galaxy. Figure 2 presents the (post-processing) Hill radius for all stars in our simulation to illustrate its variation over time.

Fig. 2. Approximate size of the Hill radius as a function of time for all the stars in the simulation since their birth until today and 100 Myr into the future. The mass of the Oort clouds is color coded (key to the right). The star HD 103095 (see also Fig. 14) is again easily distinguished by its curious orbit.

Kuiper cliff: a rather sudden drop in the density of asteroids at the far side of the Edgewordt-Kuiper belt (Edgeworth 1943; Kuiper 1951).

Oort cloud: The Oort cloud is reserved for bound objects but for which a 1 M_⊙ body passing on a hyperbolic orbit at the Hill radius distance induces a relative velocity-perturbation of δv = 10⁻⁵. Here we define the relative velocity perturbation δv ≡ |v_⊙ − v_Gal|/v_⊙ as the change in velocity of an object at apocenter in its orbit around a star due to the presence of Galactic tidal field (see Portegies Zwart & Jílková 2015). The inner edge of the Oort cloud is shown in Fig. 1 by the right-most cyan curve.

Parking zone: The region between the Kuiper cliff and the Oort cloud. The parking zone in the Solar System is relatively empty.

Rogue planet: Planet that is unbound from any star. We note that strictly speaking, the term free floating is misleading because they are still bound to the Galaxy.

Sōlus lapis (plural: sōlī lapidēs): Asteroid or comet that is not bound to any star but floating freely in the Galaxy. Asteroids to the far side of the Hill radius are generally not bound to the parent star, but have become a member of the Galaxy.

2.2. Simulation environment

The calculations in this study are performed using the Astrophysics MUltipurpose Software Environment (AMUSE for short, see Portegies Zwart & McMillan 2018), a modular language-independent framework for homogeneously interconnecting a wide variety of astrophysical simulation codes. It is built on public community codes that solve gravitational dynamics, hydrodynamics, stellar evolution, and radiative transport. The underlying codes are written in high-performance compiled languages that do not require recompilation when combined with Python scripts. The framework adopts Noah’s Arc philosophy, in which there are at least two codes that solve for the same physics (Portegies Zwart et al. 2009).

Most calculations are carried out using a combination of symplectic, direct N-body test-particle integration, and the tidal field of the Galaxy (quite similar to the analysis in Torres et al. 2020). For the former two, we use Huayno, which adopts a recursive Hamiltonian-splitting strategy (much like bridge; see Fujii et al. 2007) to generate a symplectic integrator that conserves momentum to machine precision (Pelupessy et al. 2012). For the asteroids, we employed the hold drift-kick-drift scheme in the integrator with a time-step parameter of 0.02 (this results in an integration time-step of about 10 days). We use the GPU-enabled version, which employs the Kirin GPU library (Portegies Zwart et al. 2007; Belleman et al. 2014). The direct N-body code and the Galactic potential are coupled using the hierarchical nonintrusive code-coupling strategy discussed in Portegies Zwart et al. (2020) with an interaction time-step of 100 yr.

We ignore stellar mass loss in this study. The most massive star in our sample is 1.23 M_⊙ (HD 346704), which will not leave the main sequence within ∼5.9 Gyr (for Solar metallicity). This is considerably longer than the time frame of the simulations.

2.3. Reconstructing the debris distribution near the Sun

2.3.1. The yo-yo approach

This study strives to achieve a representative distribution of interstellar asteroids in the vicinity of the Solar System. The yo-yo scheme is carried out in two steps: In the first step (yo), we integrate the known nearby stars backward through time in a representative Galactic potential. In the second step (yo), we integrate these same stars from their back-calculated location in the Galaxy forwards in time to the current epoch. In the second step, each star has a disk of planets and asteroids which are also integrated in the Galactic potential.

In the following sections, we discuss the following steps in the yo-yo procedure.

• A:

  Select the 200 nearby stars used in this study (see Sect. 2.3.2).

• B:

  Calculate their Galactic location 1 Gyr ago (see Sect. 2.3.3).

• C:

  Supply each star with planets and asteroids (see Sect. 2.3.4).

• D:

  Calculate star, planets, and asteroids forwards in time for 1.1 Gyr in the Galactic potential (see Sect. 2.3.5).

2.3.2. Selecting the 200 nearest stars

We start by selecting the 200 stars currently nearest to the Sun. The proper motions, radial velocities, positions, and masses of these stars are taken from the catalog of Torres et al. (2019a,b). The proper motion, distance, and position on the sky in this catalog are taken directly from Gaia DR2 (Gaia Collaboration 2018; Andrae et al. 2018). The radial velocities are derived by cross-matching Gaia DR2 with RAVE-DR5 (Kunder et al. 2017), GALAH DR2 (Buder et al. 2018), LAMOST DR3 (Zhao et al. 2012), APOGEE DR14 (Abolfathi et al. 2018), and XHIP (Anderson & Francis 2012). Only stars with a relative uncertainty of <20% on the parallax are selected. Regretfully, this excludes some interesting stars, such as Glise 876, which according to Dybczyński & Królikowska (2018) is a favorable candidate for the origin of ’Oumuamua.

After selecting the stars, their coordinates on the sky, distances, proper motions, and radial velocities are converted to Galactic position and velocity in Cartesian coordinates, assuming the solar position of (−8.34,  0.0,  0.027) kpc with a velocity of (11.1,  240,  7.25) km s⁻¹.

In Fig. 3 we present a projection in the Galactic plane of the relative equipotential surface of these stars. The red circle at ∼16 pc indicates the distance to the nearest 200 stars, which form the basis of our calculations.

Fig. 3. Equipotential surface of ∼795 stars within 50 pc of the Sun projected on the xy-plane in the Galactic potential. To make the potential more visible, we magnified the stellar masses by a factor 300. The Sun is indicated with the red dot in the middle of the figure. The red circle at ∼16 pc indicates the inner region for which the analysis in this paper is performed. This volume contains the 200 nearest stars (Torres et al. 2019b). The script to generate this image, plot_nearby_stars_equipotential_surface.py, and the required data (filename gaia_nearby_stars.amuse), are both available at figshare.

2.3.3. Reconstructing the stellar birth locations

As the first step in generating initial conditions, we integrate the 200 selected stars backward in time in the Galactic potential. We use the model of the Galaxy by Martínez-Barbosa et al. (2015).

This back-integration is performed for 1 Gyr. Some stars may be considerably older than that, and a few may even be younger. Therefore, these initial locations should be considered approximate and are meant to reconstruct the past several hundred million years of the stars’ evolution, rather than the full evolution since the birth of the Sun. More extended integration of the Galactic stars would eventually require us to take the dynamical evolution of the Galaxy into account. However, in our calculations, we adopted a semi-analytic dynamical model of the Galactic potential, which represents the current epoch, but which can probably not be extrapolated backward reliably beyond 1 Gyr.

After calculating all stars backward in time (for 1 Gyr) in the Galactic potential, they arrive at what we consider their initial positions and velocities. In Fig. 4 we present projected views of the Galaxy in the three fundamental Cartesian planes and a small close-up at the Sun’s location in the Galactic disk in the X − Y plane (top right). At the start of our simulations (1 Gyr ago) none of the 200 stars currently in the vicinity of the Sun were within 6 pc (see Fig. 4).

Fig. 4. Positions of the 200 stars at their birth location in the potential of the Galaxy 1.0 Gyr before the present. The solid curve shows the orbit of the Sun for the last gigayear. The colored background represents the adopted Galaxy potential: bottom left, top left, and bottom right: various Cartesian coordinates. The spiral structure is visible in the X − Y plane (lower left), whereas the other two panels give only an edge-on view of the Galaxy. The starting position of the Sun is indicated with the black studs, and its end position is market with a yellow studs. The planets and asteroids of all the nearby stars overlap at this initial snapshot because the planetary systems are smaller than the symbol-size in the panels. Top-right corner panel: a magnified view of 6 × 6 pc around the Sun. The outer circle shows the Hill radius at ∼0.65 pc, and the inner circle represents the inner edge of the Oort cloud at about ∼0.05 pc. In this initial image, none of the stars are within 6 pc of the Solar System. The color gradient in the top-right panel shows the local variation in the Galactic potential.

The back-calculations result in a unique position and velocity for each star, and by calculating them forwards again, the stars arrive at precisely their starting points; the orbital integration of the stars in the semi-analytic Galactic potential is not chaotic (but see the discussion in Sect. 4).

2.3.4. Initialization of planetary systems with asteroids

Once we have determined the position and velocity of each star 1 Gyr ago (see Fig. 4), planets and asteroids are added. The planets are generated using the oligarchic growth model (Kokubo & Ida 1998; Tremaine 2015) using the implementation of Hansen & Murray (2013). We assumed the planets to form from a disk with 1% of the stellar mass with a minimum distance of 10 au from the host star, and adopting an outer disk radius of 100 au. The oligarchic growth model then leads to three or four giant planets, dependent on the mass of the disk.

Figure 5 presents the masses of the planets around each star, and their distance to the host star. The inner-most planet is always initialized at 10 au, and the other planets follow the oligarchic growth model, resulting in circular orbits in the same plane. The outer-most planet then arrives at 60 to 76.3 au (see Fig. 5). Planet masses range from 0.35 M_Jupiter to 7.08 M_Jupiter.

Fig. 5. Initial planet mass as a function of orbital separation from the parent star. The mass of the host star is color coded (see the color bar to the right). In the two regimes, relatively low-mass (≲0.85 M⊙) stars acquire four planets according to the oligarchic growth model, whereas more massive stars only receive three planets.

Several regimes are visible in Fig. 5. These result from the Oligarchic model in which planet masses increase with distance to the host star. Stars more massive than about 0.85 M_⊙ tend to acquire three relatively massive planets, whereas the lower-mass stars receive four planets of somewhat lower mass. The cumulative distribution of the initial planet masses is presented in Fig. 6 (blue curve).

Fig. 6. Cumulative mass-distribution of the planets at the start of the simulation (blue) and at the current epoch (1 Gyr after the start of the simulations), in orange (for the bound planets) and green (for the unbound planets).

Once the planets are initialized, we add a population of asteroids as test particles in the planetary plane. Each star receives 1000 minor bodies per 1 M_⊙ of stellar mass. The pericenter distance was selected homogeneously between 5 au and 150 au with a random eccentricity of between 0.1 and 0.9. This choice prevents asteroids from settling in resonances with the giant planets, where they continue to consume computer resources which is not the focus of this study.

The conditions for asteroids are presented in Fig. 1, where each red dot represents an asteroid. This choice of relatively high eccentricities may appear strange for a debris disk, but it effectively places the asteroids in the planetary conveyor belt. In this region of parameter space, asteroids receive repeated gravitational kicks from the planets, driving them into wider and more eccentric orbits while preserving pericenter distance.

Before we start integrating each planetary system in the Galactic potential, the orbital plane is rotated to a random orientation (isotropically oriented in all fundamental angles). Each star with planets and asteroids then starts its orbit 1 Gyr ago at the calculated position in the Galaxy.

Each star with planets and asteroids is considered isolated with respect to the other stars. Since the stars are born throughout a wide area in the Galaxy (see Fig. 4) this assumption seems valid. Even today (when reaching the highest congregation) the stars are sufficiently far apart that the mutual influence of one star on the planetary system of another star is negligible. This is also illustrated in the top-right panel in Fig. 4, which shows the vicinity of the present-day Sun.

While generating the initial conditions of the planetary system, we ignore the binarity (or higher-order multiplicity) of some of the nearby stars. Several nearby stars are known to be multiple, and one of the important contributors to the local population of sōlī lapidēs, 61 Cyg A and B, is indeed a binary system. We also ignored the known planets orbiting some of the catalog stars. This poses less of a limitation to the starting conditions because the known planets tend to have much tighter orbits than we adopted here, and they hardly provide constraints on the wider orbits. We leave the improvement of these assumptions for future work.

In Fig. 2 we show an estimate of the Hill radius for each star in the simulation as a function of time. At the earliest phase, the Oort clouds are empty, which shows up as the dark blue colors for all stars. After some 100 Myr the respective Oort clouds start to be populated. As time progresses, these Oort clouds are depleted. The size of the Hill radius varies with the distance to the Galactic center. Near the Galactic center the Hill radius is small, whereas it grows when the star approaches apocenter in its orbit around the Galactic center. A smaller Hill radius, approaching perigalacticon, leads to enhanced mass loss from the Oort cloud.

To illustrate the orbits of the stars, we also present a projected view of the Galaxy with the orbits of all stars in Fig. 7. Here we adopt the same color coding for the stellar orbits to make the comparison with Fig. 2.

Fig. 7. Mass of the Oort cloud for all the stars in the simulation since their birth until today and 100 Myr into the future. The legend to the right gives the color coding in units of Earth masses. The star HD 103095 (see also Fig. 14) clearly has the most elliptical orbit.

2.3.5. Integrating the equations of motion

From the computed starting locations of the stars (see Fig. 4) and after planets and asteroids are situated around each (see Fig. 1), we integrate their positions (and those of the minor bodies) forwards in time for 1.1 Gyr. Asteroids, in this calculation feel the gravity of the star, the planets and the background potential of the Galaxy. The planets around a specific star feels their sibling planets, the parent star and the background potential. Each star feels its planets and the background galactic potential.

After 1 Gyr of integration all stars have arrived at (or at least near) the currently observed Gaia DR2 position. In Fig. 8 we show the back-and forwards calculation for two stars in our sample, Glise 143 and 61 Cygnus A. For the former, both the backward and forward calculations remain within a line width in the figure. The yo-yo approach in a time-resolved smooth potential will result in precisely the same forward and backward orbital trajectory, as can be seen in the left panel of Fig. 8, where we demonstrate this for Glise 143. However, the forward calculations do not always precisely match the backward calculations because sometimes one or more planets are ejected while integrating the stellar orbit forwards in time. The result of a planet-ejecting system can be seen in the right panel of Fig. 8, where we integrated 61 Cyg A in the yo-yo approach. The forwards and backward calculations in that case result in slightly different orbits.

Fig. 8. Examples of the Galactic orbits of Glise 143 (left panel) and the 61 Cyg binary (right). For Glise 143 and 61 Cyg A we plot the position of the star and asteroids 1.0 Gyr ago, 0.8, 0.6, 0.4, and 0.2 Gyr ago, and today. We plot the orbits of both stars 61 Cyg A (red dotted line) and 61 Cyg B (blue dotted line). Both stars start 1 Gyr ago at the big studs to the top of the figure. The planets for 61 Cyg A are plotted only for the current epoch (when the star is near the Sun, to the left). At this time, two planets still orbit the star, but the other two planets were ejected early in the calculation. By the current time shown here, the ejected planets have moved away considerably from the star. The main differences in the Galactic orbits of the two stars 61 Cyg A and 61 Cyg B (compare the red and blue dotted curves in the right-hand panel) is caused by the ejection of these two planets.

The ejection of a planet causes the center of mass of the star and planets to arrive at the anticipated Gaia DR2 location, but the forward-calculated stellar trajectory will deviate from the backward calculation due to conservation of energy and linear momentum (we discuss the consequences in Sect. 4).

61 Cyg A is one of the worst cases for which the backwards and forwards calculations do not result in identical trajectories in the Galaxy, because two planets of 1.79 M_Jupiter and 0.72 M_Jupiter were ejected ∼105 Myr after the start of the forward calculation. The massive planet is ejected first with a low relative velocity, but the lower-mass planet is ejected much more violently. This second planet is visible to the far right of Fig. 8 at a distance of about 7500 pc from the Galactic center (see also Fig. 9).

Fig. 9. Current view of the Mikly way Galaxy along the three major projections and a close-up centered around the Sun in the top-right corner (as in Fig. 4). The orange, green, and red points indicate the bound and unbound asteroids and the planets, respectively. The stars are shown as bright red points. The Sun’s orbit in the Galactic potential is presented as the thin black curve which starts at the black studs (in the north-north-west corner) and ends at the yellow dot (to the west). For an animation, see online.

3. Results

In the yo-yo approach, each star reaches its current location after first being calculated backward in time and subsequently forwards (with minor bodies). The forward calculation is performed with planets and asteroids up to the current epoch and then 100 Myr into the future. The current positions of the stars, planets, and asteroids are presented in Fig. 9. During the calculations each planetary system is treated as isolated with respect to the other planetary systems.

Our calculations start with a star with planets and a conveyor belt filled with asteroids. The distributions in semi-major axis and eccentricity of planets and asteroids are presented in Fig. 1. While integrating forwards in time, the conveyor belt is emptied and the asteroids either remain bound or become unbound. In the following two sections, we discuss both planets and asteroids bound to the parent star or unbound.

3.1. Population of bound planets and asteroids

We recognize several regions where an object can have a bound orbit around its parent star; it can be parked in the Kuiper belt or deposited in an Oort cloud. No asteroids are trapped in a resonant orbit with a giant planet.

3.1.1. Asteroids in the Kuiper belt and the Parking zone

In isolated planetary systems such as those simulated here the parking zone is expected to be deprived of asteroids. Still, 0.32 ± 0.01% of the asteroids in the simulations are deposited in the Parking zone of their respective star. About one-third of these populate the area within 100 au along the outer edge of the conveyor belt. This population is visible in Fig. 10 as the red shading between the right-most black solid curve and the right-most cyan curve. There is even an odd planet in the parking zone.

Fig. 10. Orbital parameters of planets (blue spots) and 141864 asteroids (red fuzz) at 1 Gyr after the start of the simulation. The various curves are described in Fig. 1. The size of each blue point is proportional to the mass of the planet it represents. The two black curves indicate the pericenter distances associated with these minimum and maximum initial orbital separations, indicating the limits of the conveyor belt. The minor bodies are presented as small red points. In Fig. 1 we present the initial conditions, and explain the terminology.

These parked asteroids are driven there through several scatterings with planets that have also migrated. Asteroids in the parking zone did not migrate along the conveyor belt but were scattered repeatedly by planets that themselves were scattered. As a consequence, the asteroid did not move along the conveyor belt but experienced a random walk until parked in a relatively wide (100 < a <  10⁴ au) orbit with a relatively low eccentricity (e ≲ 0.8). The right-most cyan curve in Fig. 10 indicates the outer edge of the parking zone (Portegies Zwart & Jílková 2015) and the inner edge of the Oort cloud.

3.1.2. Formation and early evolution of the Oort clouds

Once an asteroid passes the right-most cyan curve in Fig. 10 we consider it a member of the respective Oort cloud. Figure 11 presents the evolution of the fraction of asteroids in the Oort cloud of their respective host star. In 25 Myr, half the asteroids are removed from the planetary system. This high efficiency is consistent with the theoretical calculations by Peñarrubia (2019). At the age of 100 Myr the fraction of asteroids in the disk drops to ∼2%, and after 1 Gyr only 0.3% of the asteroids still orbit in a disk within ∼100 au of the parent star. All bound asteroids are eventually parked in some Kuiper belt, in the Parking zone or in the Oort cloud. The evaporation of the disk is consistent with the ∝t^−2/3 decay characteristic derived by Peñarrubia (2019, see his Fig. 8) for the Oort cloud.

Fig. 11. Evolution of the fraction of asteroids in the Oort cloud of their respective host stars (solid green curve). The dotted curves give the dispersion (one standard deviation) among stars. The blue dashed curve gives the fraction of planets that populate the Oort cloud. The ocher-colored curve gives the fraction that remains in the disk either in resonant orbits or in the parking zone (outside the perturbing influence of the planets).

Exo-Oort clouds build up on a timescale of a few 100 Myr (see Fig. 2). At the age of 10 Myr, 4% to 6% of the original disk population has been launched into the Oort cloud. By that time, most orbits are still in the orbital plane of the planets.

In Fig. 10 we present the orbital distribution of asteroids that remain bound to their respective stars. The solid green curve indicates where the orbital period, P_orb, equals the circularization-diffusion timescale by the Galactic tidal field (t_diff, Eq. (5) of Duncan et al. 1987). This line indicates the orbital parameters for which a particle tends to spend more time being perturbed by the Galactic tidal field (near apocenter) than by the Sun. The orbits of asteroids and planets to this curve’s right, are circularized and isotropized by the Galactic tidal field. Asteroids between the right-most cyan and the green curves have difficulty being circularized by the Galactic tidal field; they tend to persist in rather high eccentricity orbits. However, if an asteroid is pushed further along the conveyor belt, deep into the Oort cloud, the Galactic tidal field’s circularization process operates on a timescale shorter than the timescale on which the planets eject asteroids. Asteroids that penetrate the Oort cloud further (to even wider but bount orbits) tend to circularize more quickly. The mean eccentricity in the Oort cloud is then a function of distance to the parent star and of time.

In Fig. 12 we present the cumulative distribution of the minor body’s eccentricities in the Oort cloud, 50 Myr after the start of the simulation. Here we make the distinction between an inner and an outer Oort cloud. This distinction is made at a perturbation factor of δv = 10⁻⁴. This corresponds to a semi-major axis of about 6 × 10⁴ au (for a circular orbit), compared to the division between the outer edge of the parking zone and the inner edge of the Oort cloud, which is at δv = 10⁻⁵ and corresponds to a semi-major axis of ∼10⁴ au (right-most cyan curve in Fig. 10).

Fig. 12. Cumulative distribution of the eccentricity in the inner and outer Oort cloud at an age of 50 Myr after the start of the simulations. The additional green curve gives the final distribution of P(e) ∝ e2.75. The simulated eccentricity distribution approaches the theoretical curve at about 50 Myr.

In the first few tens of millions of years, when the Oort cloud is populated, the inner and outer Oort clouds have quite distinct distributions in eccentricity. This can be seen in the evolution of the mean eccentricity in these regions, presented in Fig. 13.

Fig. 13. Mean eccentricity of the asteroids in the parking zone, and the disk (within 100 au) as a function of time for the respective Oort clouds (inner and outer). Initially, all asteroids have relatively high eccentricity, but this rapidly drops in the first few tens of millions of years because of ejection by the planets. The Oort cloud population remains steady after about 100 Myr, whereas the more inner regions slowly converge to a similar eccentricity of ∼0.7.

3.1.3. The evolution of the isotropic Oort cloud

After about 100 Myr both eccentricity distributions approach the thermal distribution. The Parking-zone population at this time is still converging to the thermal distribution. For a more thorough discussion on the asteroid migration in the Solar System, we refer to Paper II (Portegies Zwart et al. 2021).

Once in the Oort cloud, the orbital distribution of asteroids becomes isotropic in inclination and their eccentricity drops from ≫0.9 to a mean eccentricity of ⟨e⟩∼0.7 (see Fig. 13). The Oort clouds are fully developed 200 Myr after the start of the simulation (see also Fig. 2).

After reaching a maximum, the number of asteroids in the Oort clouds starts to decay on a half-life timescale of about 800 Myr. This timescale is considerably longer than the theoretical estimates by Peñarrubia (2019), but shorter than estimates for the Sun’s Oort cloud derived by Hanse et al. (2018). Our calculations account for the variations in the smooth Galactic potential through which the star moves, and part of the discrepancy stems from the wide range in the ellipticities of the Galactic orbits of nearby stars (see Figs. 2 and 8). The orbit of the Sun in the Galactic center motivated Kaib et al. (2011) to discuss variations in the inner edge of the Oort cloud around the Sun. Our calculations (and those of Kaib et al. 2011) neglect the effect of passing stars, which tends to dramatically reduce the lifetime of the Oort cloud (Weissman 1986; Stern 1990; Matese & Whitman 1992; Hanse et al. 2018).

3.1.4. Population of bound planets

Planets experience a similar fate to the asteroids, but it is much harder to eject a massive planet than a mass-less asteroid. Therefore, most planets remain bound to the parent star, and stay within the planetary region, originally between 10 and 100 au. Figure 10 shows how the planets scatter along the curves of constant pericenter (flaring to the right) and constant apocenter (flaring to the left).

About 40% of the planets are ejected in the first 200 Myr. Eventually, after 1 Gyr, ∼24.8% of the planets remain bound to their respective stars, and ∼71.6% are ejected to become free-floating. Only ∼2% of the planets acquire a bound orbit in the Oort cloud of their host star. The orbital distribution of these is isotropic in inclination with relatively low eccentricity (⟨e⟩∼0.7, see Fig. 13). Figure 10 presents the orbital distribution of planets and asteroids that remain bound to their respective stars.

Planets in an Oort cloud are rare. In all our simulations, 10 out of ∼736 planets are parked in the Oort cloud. Only ∼5% of the stars in our simulated solar neighborhood have a planet in the Oort cloud.

3.2. Population of unbound planets and asteroids

Most asteroids and planets do not remain bound to their host star. These objects become rogue planets and asteroids.

3.2.1. Rogue asteroids and sōlī lapidēs

The vast majority of asteroids become free floating in the potential of the Galaxy. Figure 8 shows the evolution of the unbound asteroids for Glise 143 and 61 Cyg; the same is shown for HD 103095 in Fig. 14 at five (equally spaced) moments in time. At the start of the simulation, all asteroids (and planets) are bound to the star, but after 200 Myr most asteroids have already become unbound.

Fig. 14. Projected Galactic orbit for HD 103095, which lost all its asteroids after the first passage of the Galactic center. HD 103095 (Groombridge 1830) is a metal-poor star (log(g) = 4.71 ± 0.3 and [Fe/H] = −1.27 ± 0.2 dex; Karovicova et al. 2018) with a high proper motion. According to Gaia DR2 its mass is ∼0.87 M⊙ (but only 0.66 M⊙ according to Takeda et al. 2007). In our calculations this star has a rather elongated orbit in the Galactic potential. See Fig. 8 for similar orbital trajectories for Glise 143 and 61 Cyg.

Although the fraction of bound asteroids is generally small (see Fig. 11), there are always some asteroids in the Oort cloud; except for the halo star HD 103095 (Gaia DR2 4034171629042489088), which lost all its asteroids at 100 Myr after the start of the simulation and did not form an Oort cloud. The close perigalactic distance of this latter star, together with the chaotic reorganization of its planets, resulting in the ejection of the lowest-mass planet (∼1.59 M_Jupiter), made it hard to preserve any objects in an Oort cloud. Its projected orbital trajectory is depicted in Fig. 14.

Any escaping asteroid tends to remain within a certain distance of its parent star (see also Correa-Otto & Calandra 2019, who performed similar calculations for the Solar System), much in the same way as tidal tails from star clusters tend to remain along the trailing and leading orbital trajectory (Sandford et al. 2017). A similar structure is observed along the orbit of the Sagittarius dwarf (Majewski et al. 2003), Pal 5 (Odenkirchen et al. 2001) and other globular clusters (Piatti & Carballo-Bello 2020). After the first ∼200 Myr, insufficient time has past for the unbound asteroids to drift far away from the parent star, and the tidal asteroidal debris tails have a length of at most a few parsecs. After 600 Myr the debris tails have fully developed, extending over several kiloparsecs. The smearing and rotation of the debris tails is then determined by their orbital phase on the Galactic potential. For most of the orbit, the tidal tails closely follow the orbit of the parent star.

3.2.2. Rogue planets

The planetary systems in our simulations tend to be born in moderate stability, and some planets are ejected while interacting with each other or with the tidal field of the Galaxy. But this planet-ejection process is less common because of their large mass compared to that of asteroids. Some of these planet ejections can be the result of numerical errors, because these systems are highy chaotic. The largest positive Lyapunov exponent indicates a stability timescale considerably shorter than the timescale over which we integrate these systems. Although we adopted a direct N-body code (Huayno) for integrating the planetary systems, there is some secular growth of the energy error, which in combination with the coupling to the Galactic potential may lead to spurious planet ejections.

A consequence of planet–planet scattering is the relatively high probability of ejection of lower mass planets, leaving the more massive planets in perturbed (possibly inclined and eccentric) orbits (Cai et al. 2018). This is illustrated in Fig. 6, where we present the initial mass-function of planets, the final mass-function (after 1 Gyr), and the mass function of the escaped planets. The final bound planet mass is larger than the initial planet mass because of the relatively common escape of low-mass planets. Initially, half the planets have a mass of ∼2.5 M_Jupiter, whereas the mean mass of escaping planets is about ∼1 M_Jupiter.

In our calculations, the effect of the Galactic tidal field on the dynamical evolution of the planets is small, except for those stars that have a small perigalactic orbit (such as HD 103095). The strong perturbing effect of the Galactic center led Veras & Evans (2013) to argue that the Galactic-center region will have a higher proportion of rogue planets than the relative outskirts of the Solar neighborhood.

3.3. Current local Galactic distribution of asteroids

Figure 9 presents the distribution of asteroids at the current epoch 1 Gyr after the start of the simulation. It is no coincidence that at this moment, most objects are located around and near the Sun. Some structures along the Sun’s orbit in the Galactic potential and some streams perpendicular to the current Sun’s orbit are discernible. In the coming paragraphs, we further explore this population by presenting various slices in parameter space. The objective of this study is to better understand the local phase-space distribution of sōlī lapidēs.

3.3.1. Characterization of the local distribution of sōlī lapidēs

Figure 15 shows the relative cumulative contribution of the local density of interstellar objects in the vicinity of the Solar System of the 200 nearby stars. In Fig. 16 the contribution for each individual star is sorted according to its current (Gaia DR2) distance. The ten stars contributing the most to the local density of sōlī lapidēs are shown in red (using a smoothing based on a kernel density estimator, KDE), and all the others are blue. These ten stars are listed in Table 1, and contribute more than 70% to the local density. About half the stars hardly contribute to the current local density of interstellar minor bodies. Even if the star is close to the Solar System, the stream of unbound asteroids may simply miss our current location in the Galaxy, as is the case for the second to fifth nearest stars to the Sun (see Table 1), which contribute less than η Cephei or Glise 143.

Fig. 15. Relative cumulative contribution to the local density of minor bodies in the Solar neighborhood (sorted by their contribution to the local density of sōlī lapidēs). We used a KDE smoothing kernel with 1.41 pc.

Fig. 16. Individual contribution to the local density of sōlī lapidēs calculated using a kernel density estimator (KDE) with a smoothing length of 0.39 pc. The ten biggest contributors to the local density are indicated in red and identified with their popular name.

The total local density of sōlī lapidēs ∼0.027 M_Earth pc⁻³. With an estimated mass for ’Oumuamua of ∼1.2 × 10⁹ kg (Portegies Zwart et al. 2018a), the nearby stars contribute ∼1.2 × 10¹⁴ per cubic parsec to the local density of interstellar asteroids.

The local density of interstellar minor bodies was estimated by Portegies Zwart et al. (2018a) and Do et al. (2018) to be of the order of 10¹⁵ objects per pc⁻³. This density has contributions from evaporating Oort clouds or from the disruption of circumstellar debris disks. Both tend to be similar in scope in our analysis, because Oort cloud formation tends to be associated with debris-disk destruction (although less closely than the process discussed in Veras et al. 2011; Veras & Evans 2013, who assumed that the disk would be obliterated by the post-asymptotic giant-branch evolution of the parent star, rather than planet scattering).

It is hard to estimate how much stars further away contribute to the local density of sōlī lapidēs, in particular because, as demonstrated here, far-away stars may contribute considerably when the Sun happens to move through the tidal-debris tail. We then estimate that Oort cloud evaporation and debris-disk destruction contribute ≳10% to the observed population of sōlī lapidēs.

To further illustrate this, the list in Table 1 is sorted by distance to the sun (n-dist =1 for the closest star p Centauri). Some nearby stars contribute considerably to the local density, but many major contributors are further away, as far as Glise 143 at about 16 pc. We note that contribution to the local density depends somewhat on the adopted smoothing kernel for the density estimator. When we adopt a smoothing kernel of 1.0 pc instead of 0.38 pc, the star SZ UMa (the 70th nearest to the Sun) contributes almost as much to the local density of sōlī lapidēs as 61 Cyg B.

Figure 17 shows the stars that contribute most prominently to the local minor-body density (see Table 1) presented as contours, whereas the others are indicated with a color-shaded kernel-density estimator. This distribution is clumpy because of the minor bodies that escaped the gravitational pull of their parent stars with a relatively low velocity.

Fig. 17. Local projected density of asteroids in the solar neighborhood (within 30 pc) in two projections. On the X − Y plane for the upper panel is a projection in the galactic plane, whereas the bottom panel shows the X − Z projection. The kernel-density smoothing kernel used was 0.104 pc. The red shades give the projected density of asteroids around the Sun. The contours give the local density only for the asteroids that originally orbited one of the ten stars in Table 1.

The debris streams are hardly visible in Fig. 17 because they tend to extend over a longer range and the relatively high local density causes the subtleties of the extended tails to wash away in the background. However, when we present the sōlī lapidēs of the most contributing stars from Table 1 on a larger scale, the tidal debris streams become evident. This is illustrated in Fig. 18. Here the colored studs identify the minor bodies of the ten most important contributors. The Sun is at the origin because this is how the contributing stars were selected as nearest neighbors to the Sun (on which the picture was centered). The morphology and width of the streams depend on the orbit of the host star and orbital phase (see e.g., Fig. 8). A rather circular orbit around the Galactic center results in a narrow and elongated stream, as is seen for the star ϵ Eridani, whereas a more elliptical orbit, such as in the case of Glise 143, results in a broader distribution along its orbit. The interesting horseshoe shape of the tidal arms results from Glise 143 being close to apocenter in its orbit around the Galactic center: the trailing arm has the same spread as the leading arm but being located closer to the Galactic center it overtakes the leading arm.

Fig. 18. Current density of asteroids within 1 kpc of the Sun in the Galactic xy-plane and centered around the Sun (black point in the middle). The Galactic center is to the right. The ten stars that contribute most to the local density of sōlī lapidēs are indicated in color (see legend), the others as small gray dots.

3.3.2. The closest approaching objects

The majority of the Solar neighborhood stars experience their closest approach within a few hundred thousand years of the current epoch (see also Torres et al. 2019a). The asteroids originally belonging to these stars swarm in front of and behind them (as seen in Fig. 9). These tidal tails may cause asteroids to pass the Sun millions of years before or after the star has its closest approach to the Sun. We study this phenomenon to understand the extent to which other stars contribute to the local density of asteroids.

For this purpose, we analyze the same stars and minor bodies for the last 20 Myr and into the coming 20 Myr in the Galactic potential while keeping track of the closest approaches (the data are available in the data/NN directory in figshare). In this time frame, each asteroid has one closest encounter with the Sun. Sōlī lapidēs tend to arrive in families (as is also illustrated in Fig. 18). Each family originates from a single parent star. These families form from asteroids that escape the stellar gravity well, and continue in a similar orbit around the Galactic center to that of the parent star. A particular star can therefore lead to multiple close encounters with the Sun. These encounters may be spread out over the timescale it takes the Sun to cross the tidal-debris tail, and this may be several million years. During this time frame the Sun may have several encounters with different asteroids that originated from the same star. If several such interstellar asteroids were found, it would be much easier to find the source star.

In our simulations, the current dominant contributor to sōlī lapidēs seems to be Glise 752. In reality, Glise 752 may not contribute any sōlī lapidēs simply because, contrary to our assumptions, it may not have an Oort cloud. Glise 752, is a known binary system with a possible planet (Kaminski et al. 2018), but this does not make it necessarily a candidate responsible for the processes described here.

Most asteroids from the current neighboring stars tend to have their closest approach with the Solar System in a rather short period of only a few million years. Even those further away, or the high proper-motion halo stars in the sample (see e.g., Fig. 14), tend to have their close approach within a few million years of the current time. As we selected the 200 nearest stars, this is not a surprise, and reflects how the current nearest stars were selected.

The majority of stars orbit the Galactic center in the same general direction as the Sun, but a few have retrograde orbits with respect to the Sun. These latter are visible in Fig. 19 as the green points among the blue ones (and vice versa). These retrogade orbits can be identified in Fig. 20 as those with a high relative velocity with respect to that of the Sun.

Fig. 19. Moment of closest approach relative to the current epoch along the three major axes. Asteroids in the wake of the Solar System tend to have been encountered a while in the past, whereas those in front will have the closest encounter in the future.

Fig. 20. Relative velocity at closest approach vs. distance of closest approach. The individual families of objects can be identified as narrow streams to short distance. To guide the eye we overplot the actual (unsmoothed) data with small dots.

The streams of asteroids in the left and right wings in Fig. 19 are the result of families of objects belonging to specific stars that pass today’s Solar System. Asteroids released around the stars’ leading tidal tail in the potential of the Galaxy have a close encounter with the Sun at an earlier epoch, whereas asteroids that escape through the trailing Lagrangian point have their closest approach to the Sun at a later time. These streamers are even more pronounced in Fig. 21 where we present the relative velocity at the closest approach as a function of the time of the closest approach. Specific families of asteroids are now clearly visible, in particular at relatively low relative velocity.

Fig. 21. Distribution of relative velocity at the closest approach to the Sun as a function of the moment of closest approach. The color coding indicates density normalized to the maximum. The streams of objects with coherent motion in time represent families of asteroids belonging to a specific star.

In Fig. 22, we present the distance of the closest approach versus the moment of the closest approach. The concentration to the current epoch (at t = 0 Myr) along the x-axis supports the results presented in Figs. 19 and 21. This result is consistent with earlier estimates of the distance and moment of closest approaching nearby stars (Bailer-Jones 2015; Darma et al. 2019; Wysoczańska et al. 2020; Bobylev & Bajkova 2020).

Fig. 22. Distance of closest approach as a function of time since the current epoch. The color coding reflects relative density and was generated using a KDE with a dispersion of 1.41 pc.

When represented in encounter velocity relative to the Sun as a function of closest-approach distance (Fig. 20) the various families are clearly remain coherent. However, at the resolution of the simulations, the closest approach distance to a particular asteroid that originates from a particular star is merely a result of shot noise due to our relatively low resolution of 1000 asteroids per M_⊙.

3.3.3. Closest free-floating planets

In Fig. 23, we present the time of closest approach and its distance for the planets in our simulations. Whereas asteroids tend to be spread out over a much large phase-space volume, planets are confined to the current epoch within about 10 pc. Massive planets stay close to their parent stars. Although some are eventually ejected to become rogue planets, the ejection velocity is so low that they have an even higher tendency to remain than the asteroids. Asteroids show a similar over-density, as can be seen in Fig. 22, but the population is less pronounced due to the large number of sōlī lapidēs at a larger distance and the small number of asteroids that remain bound to their parent stars.

Fig. 23. Distance of closest approach as a function of time since the current epoch. The color coding shows relative density and was generated using a KDE with a disperison of 1.41 pc.

The relative proximity of free-floating planets to their parent stars is also visible in Fig. 9. Here, planets remain closer to the current Solar position in the Milky-way Galaxy than the asteroids.

4. Discussion

We simulated a hypothetical population of asteroids and planets in orbit around the 200 nearest stars to the Sun. In this calculation we ignored all known planets and binary systems among this population, but treated each star as single, orbiting in a semi-analytic potential of the Milky way Galaxy. Our calculations start 1 Gyr ago by calculating the positions and velocities of the 200 nearby stars backward with time to this epoch. From there, we surround each star with three or four giant planets in circular orbits between 10 au and 100 au, and 1000 minor bodies per M_⊙ in the conveyor belt (with pericenter between 5 au and 150 au).

4.1. Time reversibility and the uniqueness of the calculation

The integration is continued from 1 Gyr ago, through the present day, to 100 Myr into the future. This backwards-forwards yo-yo approach is time reversible and the star, after calculating into the future, arrives at precisely the same location with the same velocity as observed today.

However, in our simulations, the simulated positions and velocities of individual stars today do not always match the observations by Gaia DR2, because some of the stars lose one or more planets on the way, causing their orbits to deviate from the calculated orbit as discussed in Sect. 2.3.3. The median of the relative difference in a position perpendicular to the orbit of the star (given the current solar position, this is along the Galactic x-axis) is only 0.39 au, but along the orbit, the discrepancy is 96.7 au in the y-direction, and 451.7 au in the Galactic z-direction. Furthermore, there are a few cases for which the ejection of a planet leads to a considerable deviation of the orbit of the parent star in comparison to the adopted orbit when calculated without planets. One example is visible in Fig. 9 to the left of the location of today’s Solar System. This star, 61 Cyg A, has an orbit in the Galactic potential that forms a relatively large angle to the current Solar motion. As the dispersion in the orbital calculation (due to possible ejected planets) is largest along the orbit, this star is now located rather far from its anticipated location according to the Gaia data. Also, for other stars, the discrepancy between the center of mass and the actual stellar orbit is largest along the orbit.

4.2. The 1 Gyr time window

For the yo-yo approach, we chose a fixed time frame for the calculations of 1 Gyr. Ideally, we would have followed each star from its birth in the Galactic potential to the current time. However, this procedure would have been hindered by our lack of knowledge of the ages of all nearby stars and the long-term evolution of the potential of the Milky way Galaxy. Another complication arises from the uncertainty in the actual planet population and the presence of a debris disk. To achieve a consistent assessment of the potential contribution to the observability of the tidal asteroid tails of these stars would then require many calculations for each star (to integrate over the uncertainty in the kinematics and the orientation of the disk) and for these simulations to be run with a potential model that matches the history of the Galaxy. The first analysis would be expensive in terms of computer time, but not impossible. We simply lack the stellar ages and do not know the history of the Galaxy well enough to perform such a simulation.

The tidal tails develop in the first 100 Myr, and then slowly diffuse, causing the debris trail to grow in length. This later diffusion process is slow and develops on a gigayear timescale. For the length of the tidal-debris tails of the nearest 200 stars, the ages of the stars are not critical. This would be different if a particular star were found to be younger than about 200 Myr, in which case the star’s Oort cloud would still be in development, and the tidal tail would not yet have fully formed.

4.3. Erosion of the Oort cloud

According to calculations on the evolution of the Sun’s Oort cloud, Hanse et al. (2018) derive a half-life of to 10 Gyr. This range stems from studying Fig. 4 of Hanse et al. (2018) where erosion proceeds faster due to close encounters with other stars (Duncan et al. 2011; Torres et al. 2020).

Far-away encounters are not very important for the erosion of the Oort cloud, but the occasional close encounter induces a major degradation of the Oort-cloud population. Without those close encounters, Hanse et al. (2018) estimated the Oort cloud half-life to be of the order of to 10 Gyr, which is an order of magnitude longer than our estimate of Myr.

Part of the discrepancy can be understood from the lower stellar masses in our calculations, which measure ⟨m⟩ = 0.72 ± 0.21 M_⊙. As a result, the Hill radius of the typical nearby star is smaller than that of the Sun. It would therefore be harder for these stars to keep asteroids bound. This effect leads to an increase in the evaporation rate, but cannot explain the order-of-magnitude difference.

Another argument for our smaller Oort-cloud lifetime stems from the orbits of many of the stars, which have somewhat larger ellipticities than the Sun’s orbit around the Galactic center (see Fig. 7). This leads to a strongly varying Hill radius (see Fig. 2) and therefore stronger perturbation by the Galactic tidal field, which leads to the easier ejection of asteroids (De Biasi et al. 2015). We demonstrate that this process is responsible for the evaporation of the Oort cloud of Glise 752 in Sect. 3.3.2. On average, this process reduces the lifetime of the Oort cloud but is still insufficient to explain the discrepancy with Hanse et al. (2018).

We argue that the difference between the lifetime of the Oort cloud derived by Hanse et al. (2018) and that of our study stems mostly from how we form the Oort cloud from a population of asteroids in the conveyor belt rather than initializing them with a smoothed and virialized distribution function, as was done by these latter authors. The latter method results in a more stable Oort cloud and a lower evaporation rate.

With our derived half-life for the Oort cloud of ∼800 Myr, and current mass of 2−40 M_Earth (or 10¹⁰ to 10¹² asteroids, Weissman 1996), we derive an initial mass for the Oort cloud of 3.2−64 M_Jupiter. With an Oort cloud retention fraction of ∼0.06, this results in an initial debris disk of between 0.003 M_⊙ and 0.064 M_⊙. This value is consistent with the initial value of 0.01 M_⊙ which we adopted in our simulations. These numbers are consistent with earlier estimates of the mass of the Oort cloud (Weissman 1983; Marochnik et al. 1988; Dones et al. 2004).

Although ignored in our calculations, occasional close passages of stars turn out to be important for erosion of the Oort cloud. Therefore, these passing stars are also expected to be important for the extent and morphology of the debris streams of sōlī lapidēs. We ignored this effect in our calculations, and it is hard to anticipate, except perhaps in a statistical sense.

4.4. Interstellar visitors: ’Omuamua and Borisov

One of the objectives of this study is to obtain a better understanding of the contribution of nearby stars to possible interstellar visitors, such as 1I/’Oumuamua and 2I/Borisov. It was anticipated that these objects may have originated from the Oort clouds of other stars (Moro-Martín 2019).

In an attempt to find the origin of ’Oumuamua and Borisov, Dybczyński & Królikowska (2018) studied the stars within 60 pc (using Gaia DR2 data) but conclude that it would be difficult to find a stellar origin for the two sōlī lapidēs. The accuracy of the orbits of ’Oumuamua and Borisov and the scattering in the Galactic disk make it hard to trace back their origin beyond some 15 pc or within the last 10 Myr (Hallatt & Wiegert 2020).

Despite this earlier analysis, in which nearby stars were excluded as the source of the two observed sōlī lapidēs (Portegies Zwart et al. 2018a; Bailer-Jones et al. 2018, 2020), Dybczyński & Królikowska (2018) argued that Gl 876 might be a suitable host for ’Oumuamua. Currently, Glise 876 is at a distance of 4.68 pc, but 0.817 Myr ago it had the closest approach to the Sun at a distance of 2.24 pc. Regretfully, Glise 876 (Gaia DR2 2603090003484152064) was not part of our analysis. If indeed Glise 876 is a contributor to interstellar visiting rocks or icy bodies, we might expect multiple objects to come from the same source.

Objects such as ’Oumuamua and Borisov enter the Solar System on unbound orbits and leave again after interacting with the Sun. Both objects are represented with studs in Fig. 24, which indicates the point on the sky from which they appeared to be coming, calculated for ’Oumuamua by Bailer-Jones et al. (2018) and for Borisov by Bailer-Jones et al. (2020).

Fig. 24. Sky projection of nearby asteroids. The blue and orange points indicate the direction from which Borisov and ’Omuamua entered the Solar System from the analysis of Bailer-Jones et al. (2018, 2020), respectively. The color coding shows the relative contribution to the local density, whereas the contours provide the same information for the nine nearest stars in Table 1.

Overplotted in Fig. 24 is the distribution of all minor bodies in our calculations (red shaded region) and of the stars that contributed the most to nearby objects from Table 1 (contours). Borisov is quite far from these distributions, but ’Oumuamua is fairly close to some of the asteroids that originated from relatively nearby stars.

Multiple objects originating from the same stellar source do not necessarily come from the same direction or with a similar velocity. If we were able to argue that multiple intruders have a similar origin, it may in that case be possible to identify the family and the star from which they originated. This is also illustrated in Figs. 19 and 24, where we present the contributing streams of the stars that contribute the most to the density of asteroids.

The time ’Oumuamua has been floating around freely in space was estimated to be less than 1 Gyr (Portegies Zwart et al. 2018a; Hands et al. 2019; Hallatt & Wiegert 2020). In our calculations, tidal tails remain coherent for considerably longer. It is therefore conceivable that ’Oumuamua is still a member of the leading or trailing arm of debris around its parent star. In that case, we might expect more objects to enter the Solar System from the same Galactic direction.

To further verify that the nearby stars contribute prominently, but that far away stars may still contribute considerably, we present the evolution of the local (and maximum) density of sōlī lapidēs in Fig. 25. Today (at t = 0 Myr) the peak in asteroid density is about two orders of magnitude higher than at any other moment. Narrow peaks in time appear when the Solar System happens to move through the debris-stream of another star.

Fig. 25. Evolution of log10 density of interstellar asteroids in the Solar vicinity. Calculations were performed with a KDE smoothing kernel with a kernel bandwidth of 1.4 pc. Both curves are normalized to the same maximum of ∼0.027 MEarth pc−3 at the present epoch. We adopted a KDE smoothing kernel of 1 pc3.

Today’s high peak in the local density of sōlī lapidēs does not mean that we live at a special time; it is merely a consequence of our initial selection of the 200 nearby stars. At some points in time, we will move through the stream of a star, giving rise an enhanced encounter rate with the asteroids in its tidal streams. To what degree this results in a continuous encounter rate or in rather discrete peaks is at this point hard to assess. Based on our simulations, we cannot determine the temporal variations in this encounter rate because this depends on the extent and density structure of the individual streams. The resolution of our calculations was insufficient to determine this structure, and this would require us to take many more nearby (and further away) stars into account in our analysis.

5. Conclusions

We studied the possible evolution of the population of the 200 nearest stars at the present epoch from the Gaia DR2 catalog in the Galactic potential including a hypothetical population of initially bound planets and asteroids. This study is a continuation of the work by Torres et al. (2020), and extends the calculations therein over a longer time frame and takes the nearest 200 stars into account. The objective of these studies is to characterize the local environment of the left-overs of the planetary scattering process around these stars in order to constrain the distribution of asteroids and free-floating planets in the Sun’s vicinity. We ignore the contribution of other stars not near to the Sun. Our assumptions on the initial conditions are somewhat capricious because we do not know if these other stars have debris disks.

For the initial conditions, we assume each star to be born a billion years ago with three or four planets in circular orbits between 10 and 100 au. The planet masses range from one to seven times the mass of Jupiter. Each planetary system is further populated with 1000 point-mass asteroids per M_⊙ in the plane of the planets. The orbits have a pericenter distance of between 5 au and 150 au, and eccentricities randomly between 0.1 and 0.9. Therefore, they could be viewed as an excited population that already underwent some scattering from the planet-formation process and early dynamical evolution. Each planetary system is oriented randomly. The birth location for each star is determined from back-calculating its orbit in the Galactic potential.

After setting up the initial conditions, we numerically integrate all the stars, planets, and asteroids in the Galactic potential for 1.1 billion years ago, until the current epoch, and 100 Myr into the future. While performing these calculations, we keep track of the positions of the Sun and all other objects to find the moment of closest approach between any object with respect to the Sun.

We find that in the 1 billion years of evolution, 98% of the asteroids become unbound and turn into free-floating asteroids, or sōlī lapidēs. About 2% of the asteroids remain bound to their parent star in the form of an Oort cloud, and ∼0.31% remain bound to their parent star in stable orbit outside the range of influence of the outer-most planet, in the parking zone. A tiny fraction of asteroids, namely ∼0.1% (164), find themselves bound to their parent star outside its respective Hill radius. The fraction of planets that remain bound in the Oort cloud of their parent star is ∼2%, similar to the fraction of asteroids that remain bound. This fraction is consistent with the results of the calculations by Correa-Otto & Calandra (2019).

In the first few tens of millions of years, when the Oort clouds are filled, the inner and outer Oort clouds have quite distinct distributions in eccentricity. This changes after about 50 Myr, when the eccentricity distribution approaches the thermal distribution, as demonstrated in Fig. 12.

The distribution in eccentricities of the objects in the Oort clouds reaches a steady-state of P(e) ∝ e^2.75 after about 50 Myr. The distribution in semi-major axes is flat in log from 10⁴ au to 10⁵ au. The half-life of the Oort clouds in our calculations is 800 Myr, which is considerably shorter than found in earlier studies. This short half-life in our simulation is attributed to the correlation between semi-major axis, eccentricity, and inclination. However, some contribution also comes from the lower masses of the stars and their relatively high eccentricity in the Galactic potential.

The ten stars in the Gaia DR2 database that contribute most prominently to the probability of producing an interloper are not necessarily the closest stars. Some relatively distant stars contribute more prominently than some closer in. This probability depends on the orbit of a particular star (and its asteroids) with respect to the Sun. A star that orbits parallel to the Sun in the Galactic potential is unlikely to contribute to the local density of sōlī lapidēs whereas a crossing orbit results in a local over-density. A relatively faraway star can therefore still contribute considerably to the local asteroid density, in contrast to a more local star. These peaks in over-density of local sōlī lapidēs do not necessarily coincide with the closest approach for that particular host star, but may happen millions of years earlier or later.

The contribution to the locally observed population of sōlī lapidēs strongly depends on whether or not the Sun happens to move through the tidal asteroid tail of another star. Some of these trails are dense and narrow, and may give rise to a local enhancement in the incident rate of interstellar objects interacting with the Solar System. We derive a current local density of sōlī lapidēs from the 200 nearest stars of ∼0.027 M_Earth pc⁻³. With a typical mass comparable to that of ’Oumuamua, we then arrive at a local number density of ∼1.2 × 10¹⁴ pc⁻³, which is a factor of ten lower than earlier estimates based on the two known sōlī lapidēs. We therefore argue that the nearest few hundred stars contribute to the rate of interstellar asteroids visiting the Solar System, but they are probably not the dominant source.

We expect considerable structure in the distribution of interstellar asteroids. If observable, this structure will inform us about the Galactic potential, the orbits and ages of the nearby stars, and the internal reorganizations that their planetary systems might have endured. On these grounds, we expect considerable temporal variations in the encounter rate of sōlī lapidēs.

6. Energy consumption of this calculation

The evolution of each star (including the planets, asteroids and the Galactic tidal field) was computed in 3 to 7 days on one GPU with six cores. The total computer time spent was about 1000 h on the GPU and 6000 CPU hours. With 180 Watt h⁻¹ per GPU and 12 Watt h⁻¹ per CPU our total energy consumption for the calculations was about 250 kWh. Considering the emission intensity of 0.283 kWh kg⁻¹ (Wittmann et al. 2015) our calculations resulted in almost 900 tonnes of CO₂, quite comparable to launching a rocket into space (Portegies Zwart 2020).

Movie

Movie of Fig. 9 Access here

Acknowledgments

It is a pleasure to thank Fransisca Concha-Ramírezs, Anthony Brown, Martijn Wilhelm and Jorge Peñarrubia for discussions. Special thanks to Santiago Torres for numerious stimulating discussions, and very pleasant collaboration throughout his PhD. I thank the referee for helping considerably in improving the presentation and shaping the paper. This work was performed using resources provided by the Academic Leiden Interdisciplinary Cluster Environment (ALICE). This work would have been impossible without the following public open source packages and libraries: Python (van Rossum 1995), matplotlib (Hunter 2007), numpy (Oliphant 2006), MPI (Gropp et al. 1996; Gropp 2002) and AMUSE (Portegies Zwart et al. 2018b).

References

All Tables

All Figures

Fig. 1. Orbital parameters of planets (blue studs) and 141864 asteroids (red points) around the 200 nearest stars 1 Gyr ago. Here the right-most light-blue curve gives the approximate inner edge of the Oort cloud using a mean local stellar mass of ∼1.0 M⊙ and a mean distance according to the Sun’s Hill radius in the Galatic potential (∼0.65 pc, red vertical line). The vertical dashed and dotted curves give the smallest and largest Hill radii experienced by any star throughout the simulations. The inner light-blue curve is given for a relative perturbation of δv = 10−8, and the outer light-blue curve gives δv = 10−5. The green curve gives the distance from the host star for which the circularization timescale of the Galactic tidal field is comparable to the eccentricity pumping timescale of a test particle in the conveyor belt. The orbital configuration at the present epoch, 1 Gyr after the start of the simulation, is presented here. This figure is generated with the Python script plot_semimajor_axis_vs_eccentricity_of_planets_and_comets.py. This script, the required data, and the scripts to generate any of the figures in this manuscript including all required data are available through figshare.

In the text

	Fig. 2. Approximate size of the Hill radius as a function of time for all the stars in the simulation since their birth until today and 100 Myr into the future. The mass of the Oort clouds is color coded (key to the right). The star HD 103095 (see also Fig. 14) is again easily distinguished by its curious orbit.
In the text

Fig. 3. Equipotential surface of ∼795 stars within 50 pc of the Sun projected on the xy-plane in the Galactic potential. To make the potential more visible, we magnified the stellar masses by a factor 300. The Sun is indicated with the red dot in the middle of the figure. The red circle at ∼16 pc indicates the inner region for which the analysis in this paper is performed. This volume contains the 200 nearest stars (Torres et al. 2019b). The script to generate this image, plot_nearby_stars_equipotential_surface.py, and the required data (filename gaia_nearby_stars.amuse), are both available at figshare.

In the text

Fig. 4. Positions of the 200 stars at their birth location in the potential of the Galaxy 1.0 Gyr before the present. The solid curve shows the orbit of the Sun for the last gigayear. The colored background represents the adopted Galaxy potential: bottom left, top left, and bottom right: various Cartesian coordinates. The spiral structure is visible in the X − Y plane (lower left), whereas the other two panels give only an edge-on view of the Galaxy. The starting position of the Sun is indicated with the black studs, and its end position is market with a yellow studs. The planets and asteroids of all the nearby stars overlap at this initial snapshot because the planetary systems are smaller than the symbol-size in the panels. Top-right corner panel: a magnified view of 6 × 6 pc around the Sun. The outer circle shows the Hill radius at ∼0.65 pc, and the inner circle represents the inner edge of the Oort cloud at about ∼0.05 pc. In this initial image, none of the stars are within 6 pc of the Solar System. The color gradient in the top-right panel shows the local variation in the Galactic potential.

In the text

	Fig. 5. Initial planet mass as a function of orbital separation from the parent star. The mass of the host star is color coded (see the color bar to the right). In the two regimes, relatively low-mass (≲0.85 M⊙) stars acquire four planets according to the oligarchic growth model, whereas more massive stars only receive three planets.
In the text

	Fig. 6. Cumulative mass-distribution of the planets at the start of the simulation (blue) and at the current epoch (1 Gyr after the start of the simulations), in orange (for the bound planets) and green (for the unbound planets).
In the text

	Fig. 7. Mass of the Oort cloud for all the stars in the simulation since their birth until today and 100 Myr into the future. The legend to the right gives the color coding in units of Earth masses. The star HD 103095 (see also Fig. 14) clearly has the most elliptical orbit.
In the text

Fig. 8. Examples of the Galactic orbits of Glise 143 (left panel) and the 61 Cyg binary (right). For Glise 143 and 61 Cyg A we plot the position of the star and asteroids 1.0 Gyr ago, 0.8, 0.6, 0.4, and 0.2 Gyr ago, and today. We plot the orbits of both stars 61 Cyg A (red dotted line) and 61 Cyg B (blue dotted line). Both stars start 1 Gyr ago at the big studs to the top of the figure. The planets for 61 Cyg A are plotted only for the current epoch (when the star is near the Sun, to the left). At this time, two planets still orbit the star, but the other two planets were ejected early in the calculation. By the current time shown here, the ejected planets have moved away considerably from the star. The main differences in the Galactic orbits of the two stars 61 Cyg A and 61 Cyg B (compare the red and blue dotted curves in the right-hand panel) is caused by the ejection of these two planets.

In the text

Fig. 9. Current view of the Mikly way Galaxy along the three major projections and a close-up centered around the Sun in the top-right corner (as in Fig. 4). The orange, green, and red points indicate the bound and unbound asteroids and the planets, respectively. The stars are shown as bright red points. The Sun’s orbit in the Galactic potential is presented as the thin black curve which starts at the black studs (in the north-north-west corner) and ends at the yellow dot (to the west). For an animation, see online.

In the text

Fig. 10. Orbital parameters of planets (blue spots) and 141864 asteroids (red fuzz) at 1 Gyr after the start of the simulation. The various curves are described in Fig. 1. The size of each blue point is proportional to the mass of the planet it represents. The two black curves indicate the pericenter distances associated with these minimum and maximum initial orbital separations, indicating the limits of the conveyor belt. The minor bodies are presented as small red points. In Fig. 1 we present the initial conditions, and explain the terminology.

In the text

Fig. 11. Evolution of the fraction of asteroids in the Oort cloud of their respective host stars (solid green curve). The dotted curves give the dispersion (one standard deviation) among stars. The blue dashed curve gives the fraction of planets that populate the Oort cloud. The ocher-colored curve gives the fraction that remains in the disk either in resonant orbits or in the parking zone (outside the perturbing influence of the planets).

In the text

	Fig. 12. Cumulative distribution of the eccentricity in the inner and outer Oort cloud at an age of 50 Myr after the start of the simulations. The additional green curve gives the final distribution of P(e) ∝ e2.75. The simulated eccentricity distribution approaches the theoretical curve at about 50 Myr.
In the text

Fig. 13. Mean eccentricity of the asteroids in the parking zone, and the disk (within 100 au) as a function of time for the respective Oort clouds (inner and outer). Initially, all asteroids have relatively high eccentricity, but this rapidly drops in the first few tens of millions of years because of ejection by the planets. The Oort cloud population remains steady after about 100 Myr, whereas the more inner regions slowly converge to a similar eccentricity of ∼0.7.

In the text

Fig. 14. Projected Galactic orbit for HD 103095, which lost all its asteroids after the first passage of the Galactic center. HD 103095 (Groombridge 1830) is a metal-poor star (log(g) = 4.71 ± 0.3 and [Fe/H] = −1.27 ± 0.2 dex; Karovicova et al. 2018) with a high proper motion. According to Gaia DR2 its mass is ∼0.87 M⊙ (but only 0.66 M⊙ according to Takeda et al. 2007). In our calculations this star has a rather elongated orbit in the Galactic potential. See Fig. 8 for similar orbital trajectories for Glise 143 and 61 Cyg.

In the text

	Fig. 15. Relative cumulative contribution to the local density of minor bodies in the Solar neighborhood (sorted by their contribution to the local density of sōlī lapidēs). We used a KDE smoothing kernel with 1.41 pc.
In the text

	Fig. 16. Individual contribution to the local density of sōlī lapidēs calculated using a kernel density estimator (KDE) with a smoothing length of 0.39 pc. The ten biggest contributors to the local density are indicated in red and identified with their popular name.
In the text

Fig. 17. Local projected density of asteroids in the solar neighborhood (within 30 pc) in two projections. On the X − Y plane for the upper panel is a projection in the galactic plane, whereas the bottom panel shows the X − Z projection. The kernel-density smoothing kernel used was 0.104 pc. The red shades give the projected density of asteroids around the Sun. The contours give the local density only for the asteroids that originally orbited one of the ten stars in Table 1.

In the text

	Fig. 18. Current density of asteroids within 1 kpc of the Sun in the Galactic xy-plane and centered around the Sun (black point in the middle). The Galactic center is to the right. The ten stars that contribute most to the local density of sōlī lapidēs are indicated in color (see legend), the others as small gray dots.
In the text

	Fig. 19. Moment of closest approach relative to the current epoch along the three major axes. Asteroids in the wake of the Solar System tend to have been encountered a while in the past, whereas those in front will have the closest encounter in the future.
In the text

	Fig. 20. Relative velocity at closest approach vs. distance of closest approach. The individual families of objects can be identified as narrow streams to short distance. To guide the eye we overplot the actual (unsmoothed) data with small dots.
In the text

	Fig. 21. Distribution of relative velocity at the closest approach to the Sun as a function of the moment of closest approach. The color coding indicates density normalized to the maximum. The streams of objects with coherent motion in time represent families of asteroids belonging to a specific star.
In the text

	Fig. 22. Distance of closest approach as a function of time since the current epoch. The color coding reflects relative density and was generated using a KDE with a dispersion of 1.41 pc.
In the text

	Fig. 23. Distance of closest approach as a function of time since the current epoch. The color coding shows relative density and was generated using a KDE with a disperison of 1.41 pc.
In the text

	Fig. 24. Sky projection of nearby asteroids. The blue and orange points indicate the direction from which Borisov and ’Omuamua entered the Solar System from the analysis of Bailer-Jones et al. (2018, 2020), respectively. The color coding shows the relative contribution to the local density, whereas the contours provide the same information for the nine nearest stars in Table 1.
In the text

	Fig. 25. Evolution of log10 density of interstellar asteroids in the Solar vicinity. Calculations were performed with a KDE smoothing kernel with a kernel bandwidth of 1.4 pc. Both curves are normalized to the same maximum of ∼0.027 MEarth pc−3 at the present epoch. We adopted a KDE smoothing kernel of 1 pc3.
In the text