© The Authors 2026

1 Introduction

Oort (1950) proposed the existence of a cloud of comets that surrounds the Solar System, which was subsequently named after him. The Oort cloud is thought to contain about 10¹¹−10¹² comets (Brasser & Morbidelli 2013), with a total mass of around 2–3 M_⊕ (e.g., Francis 2005; Morbidelli 2005). Its shape is conjectured to be nearly spherical, with the cloud extending up to 0.5 pc from the Sun, limited by the Hill sphere of the Solar System (Oort 1950; Chebotarev 1965). Because these estimates of the Oort cloud properties are highly uncertain (mainly because of the lack of observations of the long-period comets), the creation, evolution, and even the existence of the Oort cloud remains a puzzle today. One of the main open questions is related to the moment of its formation.

Several works have explored the origin and evolution of the Oort cloud, mainly by performing numerical simulations. These studies follow two main ideas of formation: primordial formation and late formation. The primordial model occurs in the early stages of the Solar System when the Sun was still in its birth cluster (see, e.g., Oort 1950; Hills 1981; Heisler & Tremaine 1986; Duncan et al. 1987; Weissman 1996; Wiegert & Tremaine 1999; Dones et al. 2004; Levison et al. 2010). In this scenario, the planetesimals were scattered due to the interaction with the growing giant planets and started populating the Oort cloud. The late model assumes that the Sun has already left its birth cluster. In this alternative scenario the Oort cloud was created in the later stages of the Solar System formation possibly due to gravitational instabilities when the giant plants experienced an orbital resonance that caused the minor bodies to be ejected into almost unbound orbits (see, e.g., Levison et al. 2004; Brasser et al. 2006; Fouchard et al. 2006; Kaib & Quinn 2008; Morbidelli 2005; Duncan 2008; Brasser & Morbidelli 2013; Shannon et al. 2014; Dones et al. 2015; Nesvorný 2018; Shannon et al. 2019; Vokrouhlický et al. 2019; Wu et al. 2023, 2024).

It is now widely accepted that the Solar System was born within a stellar cluster (Portegies Zwart 2009; Adams 2010), and that the Sun remained a member of this cluster during the formation of the planets (Pichardo et al. 2012; Portegies Zwart & Jílková 2015; Martínez-Barbosa et al. 2016). Given that the typical lifetime of an open cluster is approximately 100 Myr and that stellar clusters are inherently dense environments, the likelihood of close stellar encounters during this period is significantly higher than it would be for a star evolving in isolation (Portegies Zwart & Jílková 2015). Consequently, the Oort cloud probably experienced gravitational interactions with the Sun’s stellar siblings during the early stages of its formation.

Several studies have investigated the evolution of the Solar System within dense stellar environments. For instance, Nordlander et al. (2017) examined the fate of a primordial Oort cloud across different models of the Sun’s birth cluster (low-, intermediate-, and high-mass). They considered the long-term evolution of the cloud, beginning at 100 Myr and extending through the late heavy bombardment phase (~ 500 Myr). Their findings suggest that a primordial Oort cloud is unlikely to survive in any of the modeled clusters, with the low-mass cluster providing the most favorable conditions for Solar System’s survival. This supports the hypothesis of a delayed Oort cloud formation (e.g., Brasser & Morbidelli 2013). Other authors (e.g., Eggers et al. 1997; Levison et al. 2010; Jílková et al. 2016; Pfalzner et al. 2018; Hands et al. 2019) have proposed that portions of the Oort cloud may have been captured through interactions with passing stars while the Sun remained in its birth cluster. Such stellar encounters have also been suggested as a possible origin for distant objects like Sedna. More recently, Wajer et al. (2024a,b) explored the dynamical evolution of planetesimals initially located in the Jupiter–Saturn region. Their work examined how stellar interactions in the birth cluster could lead to the formation of Sedna-like objects. Their simulations predict the existence of previously undetected classes of sednoids, including very small bodies and retrograde orbiters. The similarities and differences among these studies highlight a central challenge: the uncertainty in the initial conditions of the Sun’s birth cluster and the poorly constrained history of stellar encounters experienced by the early Solar System.

In their detailed numerical study, Portegies Zwart et al. (2021) provided a comprehensive picture of the chronology underlying the formation and evolution of the Oort cloud. Their simulations indicate that particles with semimajor axes ranging from ~100 au to several thousand au retain dynamical signatures indicative of the Sun’s origin within a dense stellar cluster of approximately a ≥ 1000 M_⊙p⁻³. Furthermore, their results suggest that most of the outer Oort cloud formed after the Solar System’s departure from its natal environment. Nevertheless, the initial formation processes were initiated by planet–disk and stellar–disk interactions occurring while the Solar System still resided within its birth cluster.

The formation and evolution of the Oort cloud offer a unique opportunity to better understand the formation of similar populations of comets in other planetary systems, which to date remain undetected. However, the detection of planetary systems hosting planets, debris disks, and exocomets (see, e.g., Barrado y Navascués et al. 1999; Lecavelier des Etangs et al. 2022; Rebollido et al. 2024; Gragera-Más et al. 2025) has opened the possibility of investigating the formation of Oort cloud-like structures and the dynamical processes that give rise to resonant bodies and interstellar comets or interstellar objects (ISOs). Motivated by these findings and by the hypothesis that the formation of the Oort cloud began during the early evolutionary stages of the Solar System while the Sun was still embedded in its birth cluster, we used the Solar System as a laboratory to investigate the processes that led to the formation of the Oort cloud-like structures in a dense environment. Additionally, we explored the broader implications for the formation of substructures such as the Kuiper Belt, the scattered disk, and interstellar comets in other planetary systems. This approach will help improve our understanding and interpretation of future observations of Oort cloud-like structures and the formation of distinct dynamical populations, including Sednoids and ISOs. These populations are likely common in planetary systems, but remain largely undetected at present.

To address this, we investigated the dynamical evolution of a debris disk under the combined influence of both single and multiple random stellar encounters, as well as gravitational perturbations from the giant planets. As a reference framework, we adopted the Solar System and constructed two analogous configurations: an Extended model and a Compact model. The Extended model is based on the present-day orbits of the giant planets, Jupiter at 5.2 au, Saturn at 9.5 au, Uranus at 19.2 au, and Neptune at 30.1 au, and features an extended debris disk spanning 40 < a 1000 au. The Compact model follows a Nice-based configuration (Gomes et al. 2005; Tsiganis et al. 2005; Morbidelli et al. 2005), where the giant planets are initially placed on more compact orbits: Jupiter at 5.5 au, Saturn at 8.1 au, Neptune at 11.5 au, and Uranus at 14.2 au. This configuration includes a more confined debris disk extending from 16 < a < 35 au.

The structure of this paper is as follows. In Sect. 2 we describe our numerical implementation. In Sect. 3 we present the results of our N-body simulations, which explore the evolution of Solar System analogues within a birth cluster environment under two modeled configurations. In Sect. 4 we discuss the effects of a dense stellar environment on the formation of transitional interstellar comets (Torres et al. 2019) and hyperbolic objects such as ’Oumuamua (Chambers et al. 2016; Meech et al. 2017). Finally, in Sect. 5 we summarize our results and present our conclusions.

2 Numerical implementation

To explore the formation of Oort cloud-like structures and interstellar comets, we examined the dynamical evolution of a debris disk subject to gravitational perturbations from planetary and stellar encounters. We developed two distinct models, the Compact and Extended configurations, each representing a Solar System analogue (SSA) composed of the four giant planets and a surrounding non-self-gravitating test particle disk. For each model, we conducted two simulation scenarios: one scenario featuring a single close stellar encounter between the SSA and another star (Sect. 2.1), represented by the Extended 1 and Compact 1 models; and another scenario involving multiple stellar encounters (Sect. 2.2), represented by the Extended N and Compact N models. Detailed initial conditions for these four simulation models are presented in Table 1.

2.1 Single stellar encounter

For the simulations and the generation of initial conditions for the Extended 1 and Compact 1 models, we employed the N-body code REBOUND (Rein & Liu 2012) with the IAS15 integrator. IAS15 is a high-order, adaptive timestep integrator chosen specifically due to its ability to achieve extremely low relative energy errors (~ 10⁻¹⁴). Despite being a nonsymplectic integrator, IAS15 demonstrates superior energy conservation compared to symplectic integrators, and its adaptive timestep capability facilitates accurate resolution of close stellar encounters and potential collisions. We carried out a total of 32 simulations, each modeling a stellar encounter with a Sun-like star. These encounters were characterized by impact parameters ranging from 70 to 400 au, inclinations (i) of 0°, 30°, 60°, and 90°, and a relative velocity of approximately 1 km/s, representative of typical conditions within open cluster environments (Binney & Tremaine 2008). The debris disk particles were represented as massless test particles (1000 per simulation), distributed randomly according to their initial semimajor axes. Each simulation ran for a total duration of 20 000 years, with the closest approach occurring at 10 000 years.

2.2 Multiple stellar encounters

To account for the effects of multiple stellar encounters, we simulated the Extended N and Compact N models using LonelyPlanets framework (Cai et al. 2015, 2017, 2018), which integrates NBODY6++GPU (Wang et al. 2015) for stellar dynamics and REBOUND (Rein & Liu 2012) for planetary dynamics within the Astrophysical Multi-purpose Software Environment (AMUSE; Portegies Zwart et al. 2009; Pelupessy et al. 2013; Portegies Zwart & McMillan 2018). The stellar environment was modeled as an open cluster using NBODY6++GPU with a Plummer sphere density distribution (Plummer 1911). The initial mass function (IMF) for the stars ranged from 0.08 to 20 M_⊙ (Kroupa 2001), comprising 2000 stars within a virial radius of 1 pc. We placed the SSA into this stellar environment using the LonelyPlanets module and evolved the combined system over a period of 100 Myr.

In the LonelyPlanets framework, SSA follow the global N-body dynamics of their host stars within the cluster. The resulting encounter statistics naturally reflect the underlying cluster demographics, which strongly favor low-mass perturbers. The Kroupa IMF yields an average stellar mass of only ⟨M⟩ ≃ 0.5–0.6 M_⊙, such that in a 2000-member cluster fewer than ~2% of stars exceed 2 M_⊙ and only ~0.1% exceed 8 M_⊙. SSAs are therefore more likely to reside in the cluster outskirts, where encounters are dominated by abundant low-mass stars and close approaches with massive perturbers remain rare (e.g., Stock et al. 2020). Moreover, massive stars evolve off the main sequence on timescales much shorter than the ~100 Myr dissolution time of a typical young cluster, leaving only a brief window for interactions with planetary systems. Mass segregation further enhances this bias by driving the most massive stars into the cluster core on timescales of a few Myr (Spurzem & Takahashi 1995; Mouri & Taniguchi 2002).

Once we set the initial conditions, we initially conducted a preliminary analysis using a grid of low-resolution simulations (200 simulations, each with 200 test particles in the debris disk) to efficiently minimize computational resources. From these preliminary results, we selected 24 systems for detailed, high-resolution simulations, each incorporating 2000 test particles in the debris disk. In every high-resolution run, the SSA was initialized at different positions within the cluster to capture the variety of local environments. The LonelyPlanets module employs a KD-tree algorithm to dynamically identify the five nearest neighbours to the SSA at each timestep. Following previous convergence tests in the LonelyPlanets framework, encounters with the five closest stars provide a good compromise between computational efficiency and accuracy in reproducing the cumulative perturbations experienced by planetary systems (e.g., Cai et al. 2019; Flammini Dotti et al. 2019; Stock et al. 2020; Veras et al. 2020).

The selected perturbers, together with the SSA, were then integrated using the IAS15 integrator (Rein & Spiegel 2014), which ensures precise modeling of their mutual gravitational interactions, and low relative energy errors over simulation timescales of approximately 100 Myr. The results for the Extended model are presented in Sect. 3.1, while those for the Compact models are presented in Sect. 3.2. To complement these cluster-based simulations, we also perform controlled single-encounter simulations with A- and B-type stars for the Compact and Extended models (Appendix A), to directly probe the impact of these rare but potentially encounters with massive stars.

3 Dynamical evolution of a debris disk in dense environments

Debris disks are the remnants of planet formation, such as the asteroid belt and Kuiper belt in the Solar System. Observational studies indicate that debris disk sizes typically span from approximately 10 au in their innermost regions to as much as 1000 au in their outer extents (e.g., Stark et al. 2009; Hughes et al. 2018; Bertini et al. 2023). The architecture and structure of these disks are shaped primarily by their dynamical histories during the early stages of their host stars’ formation, especially within the birth clusters. Given that most stars, including our Sun, originate within such stellar environments (e.g., Lada & Lada 2003; Portegies Zwart 2009; Adams 2010; Pichardo et al. 2012), these dense environments serve as ideal laboratories for investigating the formation of structures like Kuiper belt, scattered disk, the Oort cloud and the population of objects similar to Sedna. Using the numerical implementation described in Sect. 2, this section presents the results and detailed analysis of our simulations for both the Extended (Sect. 3.1) and Compact (Sect. 3.2) models.

3.1 Extended models

In this section we present the analysis and results of the Extended 1 and Extended N models (see Sect. 2 and Table 1 for detailed initial conditions). First, we analyse the dynamical evolution of an extended debris disk (ranging from 40 to 1000 au) subject to a single close stellar encounter (Sect. 3.1.1). Subsequently, we investigate the disk evolution under multiple stellar encounters in an open cluster environment (Sect. 3.1.2).

3.1.1 Extended 1

Following the methodology described in Sect. 2 and using the initial conditions outlined in Table 1, we conducted 32 simulations for the Extended 1 model. In Figs. 1 and 2 we present the most significant case, where the perturbing star passes within 300 au of the SSA.

Across all encounter angles, the majority of particles remain within the original disk, with retention rates of 73% for 0°, 76.4% for 30°, 81.2% for 60°, and 88% for 90°. These retained particles maintain near-zero eccentricities and semimajor axes between 40 and 60 au. Beyond this region, the stellar encounter induces the formation of three distinct dynamical structures, highlighted in green, blue, and gray in Fig. 1. In the green region, the particles exhibit eccentricities ranging from 0.1 to 0.9 with perihelion distances from 0 to 50 au, and are likely to become resonant bodies, analogous to trans-Neptunian objects in the Kuiper Belt. In the blue region, the particles have maximum eccentricities between 0.6 and 0.7 and perihelion distances from 50 to 100 au, with semimajor axes comparable to those of scattered disk objects in the Solar System. Finally, in the gray region, the particles acquire eccentricities from ~0.4 to 1 and perihelion distance of up to 1000 au, similar to the comets in the inner Oort cloud, and in particular to Sedna-like objects.

For a stellar encounter with angle at 0°, the inner region of the disk with particles with semimajor axes in the range 30 < a_c < 100 au gains approximately 8.5% more particles compared to the original distribution, while the inner Oort cloud region (200 < a < 10 000 au) acquires around 6.4% of the initial particles. Additionally, 20.2% of the particles become unbound (Fig. 2). When the inclination angle is set to 30°, the perturbation remains qualitatively similar to the 0°. The trans-Neptunian region (green area in Fig. 2) gains 1.7% more particles compared to the 0° encounter, while the inner Oort cloud region is less populated, with 5.2% of the initial particles acquiring semimajor axes between 200 and 10 000 au.

The dynamical effects of the close encounter decreased further for inclination angles of 60° and 90°. For a 60°, the population in the inner Oort cloud increases by 3 and 4.2% compared to the 0° and 30° cases, respectively. When the encounter angle is set to 90°, approximately 12% of the original particles acquire semimajor axes within the Oort cloud region. The transNeptunian region experiences a comparable population increase for 60° (6.8%) and 90° (5.2%) inclination angles.

The geometry of an encounter plays a crucial role in determining the fate of the debris disk, as the impact of stellar encounters is closely correlated with the inclination angle, as shown in Fig. 2 (see also Punzo et al. 2014; Pfalzner et al. 2018). Even with a relatively small impact parameter of 300 au, the inclination angle of the encounter significantly influences the resulting structure and dynamics of the particles within the system.

These differences are reflected in the distribution of particles across various regions of the system following the encounter. For an inclination of i = 0°, approximately 0.4% of the initial particles acquire semimajor axes within the planetary region (0 < a < 30 au), whereas for i = 30°, 60°, or 90°, no particles remain in this region. In the trans-Neptunian region, the number of particles increases for i = 30°, while the most effective inclination angle for populating the inner Oort cloud region is 90°. At this inclination, the stellar encounter raises the semimajor axes of 12% of the particles to the range between 200 and 10 000 au. The structure of the initial disk remains largely consistent across all inclination angles. Particles with semimajor axes between 40 and 65 au remain unperturbed, whereas those in the range of 65–200 au can attain eccentricities approaching 1.

As several studies have shown (see, e.g., Jílková et al. 2015; Cai et al. 2018; Pfalzner et al. 2018; Veras et al. 2020; Pfalzner et al. 2024), a very close encounter between a star and a planetary system can lead to the ejection and potential capture of a significant fraction of particles from a debris disk. Such encounters also form a population of highly elliptical objects, known as transitional interstellar objects (TIOs) (Torres et al. 2019), which may eventually evolve into ISOs due to the influence of galactic tides and subsequent stellar interactions.

Our simulations indicate that the formation of TIOs and ISOs is strongly dependent on the inclination angle of the stellar encounter, as we can see in Fig. 2. The production of ISOs changes with the encounter angle: for 0°, approximately 20.2%, while for 30° and 60°, the fractions decrease to 18.4 and 9.4%, respectively. Interestingly, for an encounter angle of 90°, none of the particles became unbound from the system despite the close encounter.

Fig. 1 Extended 1 model. Semimajor axis as a function of eccentricity and particle count after an encounter with a 1 M⊙ star at 300 au. The color bar indicates the perihelion distance of the particles. Each panel corresponds to a different encounter inclination angle (0°, 30°, 60°, and 90°).

Fig. 2 Extended 1 model. Cumulative distribution of the final energy of particles in the disk. The colored lines correspond to different inclination angles of the encounter (0°, 30°, 60°, and 90°), while the black dotted line represents the initial particle distribution as defined in Table 1. Negative energy values indicate particles that remain bound to the system after the encounter, whereas positive values correspond to ISOs. The Kolmogorov-Smirnov (KS) probabilities relative to the initial distribution for the different angles are 0.90% for 0°, 0.19% for 30°, 0.07% for 60°, and 0.05% for 90°.

3.1.2 Extended N

Following the methodology described in Sect. 2.1 and using the initial conditions outlined in Table 1 for the Extended N model, we conducted a set of 200 simulations. Each simulation accounts for the effect of the planets and the five closest stars in the cluster, as described in Sect. 2.2, on the SSA’s debris disk at any given time over a total integration period of 100 Myr.

Figure 3 shows the final distribution of the semimajor axis as a function of eccentricity and inclination for all particles across the 200 simulated systems. The figure highlights the emergence of structures reminiscent of the outer Solar System. Particles within the range of 30 to 40 au become trapped in mean-motion resonances with Neptune (2: 3, 3: 5, 4:7, 1:2, and 2:5), while particles between approximately 40 and 100 au attain eccentricities of up to 0.4 and inclinations reaching 20°. These distributions closely resemble the orbital characteristics of objects in the Kuiper belt and the scattered disk of the Solar System.

The inner region of the disk in all simulations experiences sufficient perturbations to develop substructures. In the outer disk, approximately 7.5% of the particles are ejected, significantly contributing to the population of the inner Oort cloud (1000 < a < 20 000 au). About 18% of the particles acquire highly eccentric orbits (e > 0.8); a fraction of these particles develop planet-crossing perihelia, potentially undergoing additional scattering events that place them on wider orbits with even higher eccentricities. Only a very small fraction (~ 0.02%) reaches the outermost regions consistent with the current outer Ort cloud (a > 20 000 au). Furthermore, a substantial fraction of particles (~ 36%) is ejected from their systems due to cumulative gravitational interactions with passing stars, thereby contributing significantly to the population of ISOs within the cluster environment.

We now examine in detail a selection of individual systems as shown in Figs. 4, 5, and 6. From the set of 200 simulations, we selected five representative systems (128, 14, 100, 86, and 157) that experienced a range of stellar encounter intensities, characterized as weak (systems 128 and 14), moderate (system 100), and strong (systems 86 and 157). In Fig. 4, we present the evolution of the semimajor axis as a function of time, alongside the stellar encounter distances, highlighting the closest stellar approaches (top panels of Fig. 4) for each of the five selected systems. In Fig. 5, we show the final distributions after 100 Myr of the semimajor axis, eccentricity, inclination, perihelion, and aphelion distances, along with their respective particle energy distributions. Finally, in Fig. 6 we present the resulting phasespace distributions of the particles, highlighting the structures formed within the disk due to stellar encounters.

The most stable system (weakly perturbed) is no. 128. All the particles remain in the disk after multiple encounters. The five closest encounters over 100 Myr range from 1000 to 5000 au (first panel, Fig. 4) are not close enough to inflict a considerable perturbation. However, the disk is heated somewhat in its outer parts, with a small fraction of particles reaching semimajor axes up to 2000 au. In system 14 (second row in Fig. 5), approximately 99.5% of the initial particle population remains bound; however, an early close stellar encounter (second panel, Fig. 4) causes 2.1% of particles to acquire eccentricities between 0.4 and 1, and 7.1% to reach semimajor axes greater than 1000 au. System 100 is moderately perturbed (third row, Fig. 5), roughly 71% of particles retain semimajor axes between 40 and 1000 au, while 9.6% of particles migrate outward, achieving semimajor axes between 1000 and 5000 au and eccentricities in the range 0.1–1. Additionally, most particles (78.4%) maintain orbital inclinations below 50°, whereas the remaining inclinations span from 50° to 150°. Approximately 17% of the particles are ejected from this disk, primarily due to a very close encounter (~ 300 au) with a low-mass star (~ 0.2 M_⊙).

The most strongly perturbed systems are systems 86 (fourth row in Fig. 5) and 157 (fifth row in Fig. 5). System 86 experienced two very close stellar encounters: first with a red dwarf star (0.5 M_⊙) at a distance of 488 au, followed later by an encounter with a brown dwarf (0.08 M_⊙) at 899 au. These encounters resulted in the ejection of approximately 13.4% of the initial particles, subsequently forming three distinct particle populations. The first population (green shaded region, fourth row of Fig. 5) comprises particles near 200 au with eccentricities ranging from 0.2 to 1. The second population (blue shaded region, fourth row of Fig. 5) occupies semimajor axes between 300 and 500 au, with orbital inclinations predominantly between 0° and 20°. About 23% of these particles attained eccentricities between 0.4 and 1, while the remainder had eccentricities spanning from 0.2 to 1. Finally, the third population (gray shaded region, fourth row of Fig. 5) consists of approximately 2.1% of the initial particles distributed across semimajor axes ranging from ~500 to 6000 au, eccentricities from 0.2 to 1, and orbital inclinations up to 80°.

Among all our simulations, system 157 is the most dynamically perturbed and develops orbital distributions closely resembling those observed for trans-Neptunian objects and long-period comets in the Solar System. Two very close stellar encounters at impact parameters of 336 and 507 au (Fig. 4), involving stars with masses of 0.5 M_⊙ and 0.1 M_⊙, respectively, result in the ejection of the system of approximately 57.5% of the initial particle population. Similar to the previously discussed systems, three primary populations and multiple secondary branches are subsequently formed due to weaker stellar encounters (fifth panel of Fig. 5). The first population (green shaded region, fifth row of Fig. 5) has semimajor axes between 40 and 100 au, eccentricities ranging from 0 to 0.99, and orbital inclinations from 5° to 50°. Notably, particles with semimajor axes of approximately 40–50 au enter mean-motion resonances with Neptune. The second population (blue shaded region, fifth row of Fig. 5) occupies semimajor axes between 100 and 300 au, eccentricities from 0.3 to 0.99, and inclinations up to 20°. The third and final population (gray shaded region, fifth row of Fig. 5) spans a broader range of orbital inclinations, from 20° to 180°, with eccentricities ranging between 0 and 0.99. Approximately 2.2% of the particles reach semimajor axes between 1000 and 4000 au, populating the inner Oort cloud region. Despite the strong encounters experienced by this system, only a small fraction of particles reach the outer Oort cloud region (a_c > 10 000 au), primarily because the closest stellar encounters occurred within the inner disk region, resulting in a substantial ejection of particles.

In all five selected systems, a tail emerges in the semimajor axis–eccentricity parameter space due to a combination of strong stellar encounters followed by weaker perturbations. In certain cases (systems 100, 86, and 157), multiple distinct branches appear, primarily resulting from subsequent weaker stellar encounters. Additionally, three characteristic particle populations are clearly identifiable (highlighted in green, blue, and gray shaded regions in the first column of Fig. 5). The repeated stellar encounters trigger particle ejections extending to the boundary of the present-day inner Oort cloud region (~ 20 000 au). Subsequently, perturbations induced by the Galactic tidal field will circularize the orbits of particles with large semimajor axes (see Brasser et al. 2012), eventually transforming these particles into stable Oort cloud objects. Conversely, particles maintaining smaller semimajor axes, where the Galactic tidal field is ineffective, remain on (highly) eccentric orbits, rendering their orbits dynamically unstable over secular timescales. Thus, while the observed tail structure is transient, it will ultimately evolve into a stable population of objects sharing similar orbital characteristics, analogous to the trans-Neptunian objects and Sednoids observed within our Solar System. The overall dynamical architecture of each system strongly depends on the geometry, stellar mass, and velocity of each encounter, as previously discussed in Sect. 3 and illustrated in Fig. 6.

Interestingly, after 100 Myr, the planets remained in stable orbits in all five of the analysed systems. This result indicates that, even within dense stellar cluster environments, planetary systems may endure multiple close stellar encounters without significant orbital modifications of their planets.

In summary, the simulations presented in this section illustrate different scenarios of how planetary systems can be dynamically perturbed during their early evolutionary stages. A key condition for long-term stability is that the disk of planetesimals must initially reside sufficiently far from the planets. If the planetesimal disk is too compact or located too close to the planetary region, planet–disk interactions can trigger planetary orbital instabilities, dynamically heating the disk and resulting in the ejection of a large fraction of particles to interstellar distances, thus contributing significantly to the population of interstellar comets. In the scenarios explored here, the formation of Oort cloud-like structures could have been initiated solely by stellar encounters without invoking planetary orbital instabilities. This scenario assumes that the planets already occupied their current positions while the planetary system was still embedded within its natal stellar cluster and that the protoplanetary gas disk had already dissipated (see also Brasser et al. 2012). Subsequent dynamical evolution, driven by perturbations from the Galactic tidal field and occasional additional stellar encounters, will facilitate the gradual transfer of particles from the inner Oort cloud region to populate the outer Oort cloud (Portegies Zwart et al. 2021).

Fig. 3 Extended N model. Semimajor axis as a function of eccentricity (top panel) and orbital inclination (bottom panel) for all the particles in the 200 simulated systems. The gray areas represent the different regions of the Solar System: Kuiper belt (KB), scattered disk (SD), Sednoids, inner and outer Oort cloud (OC); the dashed lines show the 2:3, 1:2, and 2:5 resonances with Neptune. The integration time is set to 100 Myr. An animation is available online.

Fig. 4 Extended N model. Orbital evolution of particles in the disk over 100 Myr. The bottom panels in each plot depict the semimajor axis as a function of time, while the top panels show the distance of the perturber. The panels represent systems 128, 14, 100, 86, and 157 (as labeled). The colored lines correspond to individual particles within each system.

Fig. 5 Extended N model. Orbital elements of particles in the disk after 100 Myr and multiple stellar encounters. First column: semimajor axis as a function of eccentricity, with particles color-coded by their perihelion distances. The green, blue, and gray shaded regions highlight different populations formed due to stellar encounters. Second column: perihelion as a function of orbital inclination, with dots color-coded by aphelion distance. The gray shaded areas represent the distinct regions of the Solar System: Kuiper Belt (KB), scattered disk (SD), Sednoids, and the inner Oort Cloud (OC). Third column: distribution of orbital energy for the particles. The red histograms correspond to the initial energy distribution, and the blue curve represents the final energy distribution. The blue curves with positive values indicate interstellar comets. Each row corresponds to systems numbered 128, 14, 100, 86, and 157, respectively.

Fig. 6 Extended N model. Phase-space distribution (X vs. Y) of particles in the disk after 100 Myr. The color scale indicates the particles’ positions along the vertical (Z) axis. The panels (from top to bottom) show systems 128, 14, 100, 86, and 157.

3.2 Compact models

For the construction of the Compact 1 and Compact N models, we followed the classical Nice model (Gomes et al. 2005; Tsiganis et al. 2005; Morbidelli et al. 2005) for the early Solar System. Although newer versions of the model have been proposed (Nesvorný & Morbidelli 2012; Morbidelli & Nesvorny 2020), Morbidelli et al. (2007) noted that the orbital evolution of the planets remains largely consistent with the original version. The classical Nice model, therefore, serves our purpose by providing a framework to understand the dynamical instability of the giant planets and the subsequent scattering of planetesimals, leading to the formation of Oort cloud-like structures in a compact planetary configuration. In Sect. 3.2.1 we present our results for a single encounter, while in Sect. 3.2.2, we analyse the outcomes of multiple stellar encounters.

3.2.1 Compact 1

Following the methodology described in Sect. 2 and adopting the initial conditions listed in Table 1, we performed 32 simulations for the Compact 1 model, varying the inclination of the encounter (0°, 30°, 60°, and 90°). We find that the minimum pericenter distance at which the SSA becomes significantly perturbed and begins to produce ISOs is q_⋆ ≃ 75 au. Figure 7 shows that the efficiency of disruption is strongly geometry dependent. Coplanar, prograde flybys maximize the energy transfer to disk particles, whereas inclined or polar encounters are comparatively inefficient.

In the coplanar case (i = 0°), 2.4% of particles are ejected from the system, and the post-encounter disk develops three distinct dynamical families. A small fraction forms a scattered population (~ 6%), characterized by low perihelia and high eccentricities (e ≃ 0.4–1), with Tisserand parameters relative to Uranus below the scattering boundary (T_U ≲ 2.5). These values place the particles deep within the Uranus-coupled regime, where high relative velocities lead to strong, chaotic encounters and efficient diffusion in orbital energy and angular momentum (Carusi et al. 1987; Levison & Duncan 1994). A second group, the scattered tail (~ 5.5%), occupies a narrow band near the Uranus-scattering boundary (2.5 ≲ T_U ≲ 3). These objects were shaped by earlier scattering events and remain marginally coupled to Uranus. Although close to the nominal decoupling limit, their relative velocity is sufficiently low to allow continued, though weaker, gravitational encounters (Valsecchi & Manara 1997), leading to moderate excitation in eccentricity (e ≃ 0.2–0.5) and inclination. The majority of particles (~ 88.5%) form a lower tail with T_U ≳ 3, corresponding to orbits dynamically decoupled from close-encounter scattering with Uranus, though still subject to secular and resonant perturbations. These objects remain only weakly perturbed from their initial configuration. Overall, the wide range of Tisserand parameters (0.98 < T_U < 8.3) reflects the coexistence of strongly scattered, marginally coupled, and dynamically decoupled populations produced by the stellar flyby.

The Tisserand parameter, T_U, serves as a quasi-invariant diagnostic of the restricted three-body problem. It constrains how a particle’s orbital elements, semimajor axis, eccentricity, and inclination evolve during gravitational encounters. Lower values of T_U correspond to more chaotic, dynamically heated trajectories (e.g., Murray & Dermott 2000). At higher inclinations, the morphology of the disk changes, but the same three families persist. At i = 60°, no particles are ejected, the scattered population nearly disappears (≃1.7%), and the scattered tail reaches ≃6.3%, leaving about ≃92% of the disk in the lower tail. For i = 90°, the scattered population vanishes entirely, the scattered tail decreases to ≃3%, and the majority of the particles in the disk (≈97%) survives in a dynamically cold state.

These results confirm that the three-branch tail seen in Fig. 7 is a hallmark of coplanar or moderately inclined encounters, where gravitational focusing is strongest. In contrast, polar flybys leave the primordial disk largely intact, producing only weak heating without significant ejection.

The cumulative energy distributions (Fig. 8) confirm these trends. Only the coplanar case (0°) shows a pronounced unbound tail (2.4%), with a smaller contribution at 30° (0.5%). At 60° and 90°, the CDFs overlap almost perfectly with the initial distribution, indicating that such geometries do not efficiently generate interstellar comets. The step-like features in the CDF mirror the three families identified in Fig. 7: the scattered population produces the most energetic tail, the scattered tail contributes a secondary rise at intermediate energies, and the lower tail represents the bulk of the disk that remains weakly perturbed and bound. Together, Figs. 7 and 8 show that close solar-mass flybys at q_⋆ ~ 75 au can only weakly populate an Oort-like reservoir, and that the efficiency declines sharply with increasing inclination for prograde geometries.

In contrast, retrograde encounters are expected to be less effective in disrupting the disk. Because the perturber and disk particles move in opposite directions, their higher relative velocity shortens the interaction time and suppresses gravitational focusing, thereby reducing angular-momentum exchange and limiting both ejection efficiency and dynamical heating of the disk (e.g., Winter et al. 2018).

Fig. 7 Compact 1 model. Semimajor axis as a function of eccentricity after an encounter with a 1 M⊙ star at q⋆ ≃ 75 au. The color bar indicates the perihelion distance of the particles. Each panel corresponds to a different inclination angle of the encounter (0°, 30°, 60°, and 90°).

Fig. 8 Compact 1 model. Cumulative distribution of the final orbital energy of disk particles. The colored lines correspond to different inclination angles of the encounter (0°, 30°, 60°, and 90°), while the black dashed line represents the initial distribution. Negative values (−1/a < 0) indicate bound orbits, while positive values correspond to ejected ISOs.

Fig. 9 Compact N model. Semimajor axis as a function of eccentricity (top panel) and orbital inclination (bottom panel) for all the particles in the 200 simulated systems. The gray areas represent the different regions of the Solar System: Kuiper belt (KB), scattered disk (SD), and the Sednoids. The dashed lines show the mean motion resonances with Neptune. The integration time is set to 100 Myr. An animation is available online.

3.2.2 Compact N

As described in Sect. 2, we performed a set of 200 simulations for the Compact N model, set up analogously to the Extended N models. Each system evolved for 100 Myr, accounting for multiple stellar encounters throughout the evolution. Figure 9 shows the final orbital distribution of all particles. Particles with semimajor axes between 15 and 20 au are found in mean-motion resonance with Uranus, resembling Kuiper Belt objects in the Solar System. The outer part of the disk experiences only mild perturbations due to the stellar encounters. Approximately 0.8% of particles are scattered into the Kuiper Belt region, 0.45% into the scattered disk, and 0.40% into the Sednoid region. Notably, none of the particles reach the Oort Cloud region. Overall, the strongest perturbations experienced by the particles arise from interactions with the giant planets, particularly Uranus. After 100 Myr, the orbital architecture of the giant planets remains stable. No planetary migration is observed in any of the simulations, primarily due to the distant nature of the encounters. Around 12.5% of the particles acquire eccentricities between 0.1 and 0.7, while approximately 7% are ejected from the system, forming interstellar comets.

In Sect. 3.1.2 we studied the five most interesting cases in our Extended N simulations. In this section, for comparison, we analyse the same systems, i.e., numbers 128, 14, 100, 86, and 157. In Figs. 10, 11, and 12 we show the evolution of the orbital elements of the disk particles, the final orbital elements and the energy distribution, and the position of the particles after 100 Myr, respectively. System 128 (first row in Fig. 11) has an early close encounter (1000 au), which in combination with the interaction between the giant planets and the disk and the following encounters leads to the ejection of 6.5% of the initial particles. The remaining particles concentrate in the original disk (95%) with eccentricities near zero, while the rest reach the scattered disk region with eccentricities between 0.2 and 0.6. System 14 (second row in Fig. 11) does not suffer an important change in the orbital parameters and remains mostly stable during 100 Myr. A small fraction (1.5%) of particles reach high eccentricities (0.1 < e < 0.6), and approximately 6% are ejected mainly because of the interaction with Uranus. The particles between 16 and 20 au enter in mean-motion resonance with Uranus. Systems 100 (third row in Fig. 11) and 86 (fourth row in Fig. 11) retain 94.5% of the initial particles.

The most significant perturbation is due to Uranus, which perturbs all the particles between 16 and 20 au (as in systems 14 and 128). Such particles become resonant bodies with Uranus. In systems 100 and 86, about 5.5% of the particles end up on hyperbolic orbits, while the two systems have a similar disk structure after 100 Myr (Fig. 12). However, system 100 suffered only one close encounter with a very low-mass star (0.3 M_⊙) at 273 au at an early stage (~ 20 Myr,), while system 86 faced two close encounters, the first one with a red dwarf star with a mass of 0.5 M_⊙ at 488 au and later on with a brown dwarf (0.08 M_⊙) at 899 au. The most perturbed system is 157. Approximately 7% of the particles are ejected from the system, attaining hyperbolic orbits. This system faced two very close encounters at 336 au and 507 au, with masses of 0.5 M_⊙ and 0.1 M_⊙, respectively. Both encounters took place at an early stage (10 and 27 Myr). Even though system 157 faced two close encounters, the particles with semimajor axes between 20 and 36 au do not suffer an important change.

The dynamical evolution of the particles in a compact disk depends primarily on the evolution of the giant planets. However, close encounters define their further evolution depending on the mass of the perturber. In the five systems presented here, we examined different scenarios. For stellar encounters, we find that in all the cases, a small fraction of the particles were ejected from the systems (~ 5–7%) in hyperbolic objects. The number of objects ejected not only depends on the perturbation due to the planet closest to the disk, but also depends on the proximity of the encounter. This suggests that the SSA faced at least one encounter in its early evolution, which partly determined the evolution of the inner regions. It is essential to highlight that in all our experiments, the inner SSA, particularly the planets, do not get perturbed even when a close encounter occurs. The scenario presented here is a simplification of the classic Nice model, which, however, gives us an idea of the importance of the giant planets in forming and shaping a Kuiper belt-like structure in planetary systems.

Fig. 10 Compact N model. Orbital evolution of particles in the disk over 100 Myr. The bottom panels in each plot depict the semimajor axis as a function of time, while the top panels show the distance of the perturber. The panels represent systems 128, 14, 100, 86, and 157 (as labeled). The colored lines correspond to the individual particles within each system.

Fig. 11 Compact N model. Orbital elements of particles in the disk after 100 Myr and multiple stellar encounters. First column: semimajor axis as a function of eccentricity. The particles are color-coded by their perihelion distances. The green, blue, and gray shaded regions highlight different populations formed due to stellar encounters. Second column: perihelion as a function of orbital inclination. The dots are color-coded by aphelion distance. The gray shaded areas represent distinct regions of the Solar System: Kuiper Belt (KB) and scattered disk (SD). Third column: distribution of orbital energy for the particles. The red histograms correspond to the initial energy distribution, and the blue curve represents the final energy distribution. The blue curves with positive values indicate interstellar comets. The rows (from top to bottom) correspond to systems 128, 14, 100, 86, and 157.

Fig. 12 Compact N model. Phase-space distribution (X vs. Y) of particles in the disk after 100 Myr. The color scale indicates the particles’ positions along the vertical (Z) axis. The panels (from top to bottom) show systems 128, 14, 100, 86, and 157.

3.3 Extended versus Compact models

In Sects. 3.1 and 3.2 we presented our results for the numerical simulations for the Extended and Compact models, respectively. First, we focused on the effects of a single but close encounter (Sects. 3.1.1 and 3.2.1) and then studied the effect of multiple encounters (Sects. 3.1.2 and 3.2.2) in a star cluster environment, over a time span of 100 Myr. For the Extended 1 model, we found that a close encounter with the present-day planet configuration of the SSA and a disk extending up to 200 au can produce different structures and populations in the disk, but the planets remain in stable orbits. On the other hand, for the Compact 1 model, a close encounter triggers the migration and excitation of the planets and the debris disk. The original disk gets perturbed, but most of the particles remain in the system. If we compare Figs. 2 and 8, we see that the production of unbound particles is higher for the extended case (5.6% of the original disk) than for the compact model (1.4%). It is important to highlight that even when a perturbing star comes close (~ 300 au) to a planetary system, the evolution of the planetesimal disk will be determined by its size. If the disk is compact, the perturbations due to the planets will dominate over the effect of the stellar encounter. If the disk is large, the planets will not suffer any orbital variation, while the particles in the disk can be heavily perturbed.

When considering multiple stellar encounters (Figs. 3 and 9) the effect on the disk for the case of the Extended N model is similar to the case of one encounter. However, the multiple effects of the passing stars create several populations of particles. These structures can be associated with different families in the Solar System, particularly the Sednoids region (Jílková et al. 2016; Pfalzner et al. 2018). This region (100 < a < 1000 au; Fig. 3) is populated with large numbers of particles with eccentricities ranging from 0.2 to 0.9 with orbital inclinations from 25 to ~150 degrees. When we compare this to the Compact N model, we can see that the disk is barely perturbed due to the passing stars. The major perturbations are because of the interaction with the giant planets, which produces at least six resonances with the particles with semimajor axes between 16 and 20 au (Fig. 9).

Overall, our simulations indicate that stellar encounters, whether single or multiple, distant or close, represent the most efficient mechanism for populating the inner Oort cloud-like region in planetary systems. Such perturbations result in a significant and rapid outward transport of particles from the planetesimal disk, typically on timescales of a few megayears. Conversely, giant planets primarily influence and sculpt the inner small SSA regions over comparatively longer periods, on the order of tens of megayears. Our results further suggest that Kuiper belt-like structures originate predominantly from dynamical instabilities and gravitational interactions with giant planets, in agreement with various realizations of the Nice model (e.g., Morbidelli et al. 2007). Extending these findings, the overall architecture and orbital configuration of planetary systems are ultimately determined by their birth environment, initial disk size, dynamical evolution, and planetary configurations. Given that planet formation is considered a common process around stars, structures analogous to the Oort cloud and Kuiper belt should be universal (see also Hands et al. 2019). These features are shaped primarily by stellar birth environments, which consequently dictate their long-term stability and survivability. Furthermore, this scenario implies that the ejection of interstellar comets is likely a widespread phenomenon, resulting in interstellar space being populated by numerous objects expelled from their parent planetary systems.

Both the Extended and Compact models show that the disk’s vertical structure is governed primarily by the encounter geometry. Coplanar and moderately inclined flybys (i ≲ 30°) efficiently excite eccentricities and semimajor axes while preserving a flattened configuration with inclinations below ~40°. In contrast, high-inclination and polar encounters (i ≳ 60°) inject vertical angular momentum, broadening the inclination distribution and producing thicker, more isotropic remnants. These results explain how early stellar flybys can create dynamically excited yet geometrically thin disks, while repeated encounters and planet–disk interactions progressively isotropize the outer regions and leave the inner disk dynamically cold (see Figs. 3, 5, 9, and 11). This mechanism supports the idea that the inner Oort cloud might have a disk-like origin (e.g., Fouchard et al. 2018; Fouchard et al. 2020, 2023), while the outer Oort cloud will become isotropic due to repeated Galactic tides and stellar perturbations (Portegies Zwart et al. 2021).

4 Interstellar comets

The discovery of the first interstellar comet, 1I/’Oumuamua (e.g., Chambers et al. 2016; Meech et al. 2017), and the subsequent detections of 2I/Borisov (Borisov 2019; Guzik et al. 2019) and 3I/Atlas (Seligman et al. 2025), have opened a new era in the exploration and understanding of the role of comets in the evolution of planetary systems. The origin of such ISOs, whether asteroidal or cometary, remains an active area of research, with various formation scenarios and dynamical evolution pathways proposed (see, e.g., Portegies Zwart et al. 2018; Torres et al. 2019; Hands et al. 2019; ’Oumuamua ISSI Team et al. 2019; Pfalzner et al. 2021; Jewitt 2024; Pérez-Couto et al. 2025). Among these scenarios, the most extensively studied involve gravitational perturbations caused by massive bodies, such as nearby planets, stellar encounters, or combinations thereof. For example, Hands et al. (2019) demonstrated that planetesimals can escape from stellar clusters with terminal velocities (v_∞) on the order of a few km s⁻¹, implying the existence of numerous objects similar to ’Oumuamua, traversing interstellar space after being ejected from their birth clusters. Furthermore, Torres et al. (2019) showed that comets located in the outermost regions of planetary systems may become TIOs, objects dynamically perturbed onto interstellar trajectories by secular gravitational interactions with distant passing stars. Such objects can pollute the solar neighborhood and, in some cases, potentially collide with nearby planetary systems (Torres, S. et al. 2023). In this section we analyse the unbound particles produced by planet–disk and stellar–disk interactions in our simulations, focusing on both the Extended N and Compact N models described in Sects. 3.1 and 3.2, respectively.

Our simulations indicate that the hyperbolic ejection velocities (v_∞) of particles escaping from the simulated SSA for the Extended N and Compact N models typically range between 1 and 3 km/s (Fig. 13). Particles in the Extended N model acquire greater kinetic energy due to strong gravitational perturbations caused by passing stars, achieving ejection velocities ranging from ~0.5 to approximately 40 km/s (inset, top panel of Fig. 13). Approximately 36% of the total particle population across all 200 simulated systems become interstellar comets (Fig. 13). Conversely, in the Compact N model, planet–disk interactions exert weaker perturbations compared to those caused by stellar encounters, resulting in roughly 7% of particles being ejected with velocities between 0.5 and 40 km/s (inset, bottom panel of Fig. 13). Overall, the ejection of asteroid- and comet-like objects into interstellar space within a cluster environment is predominantly driven by gravitational interactions with passing stars, especially when encounter angles are within 0° to 30° (Fig. 2). If planetary systems near the disk undergo instability, inducing significant planetary migration via planet–disk interactions, the efficiency of ISO production due to planetary perturbations may surpass that of stellar encounters (Raymond et al. 2010). However, in systems with stable planetary orbits, only a small fraction (~ 7%) of the particles is typically ejected into interstellar space.

The subsequent evolution of interstellar comets will primarily be governed by dynamical interactions with other stars within their birth cluster, and those objects that ultimately escape the cluster will further experience perturbations induced by the Galactic tidal field. These long-term gravitational effects will increase their hyperbolic excess velocities from a few km/s up to tens of km/s, similar to the observed velocity of ’Oumuamua. Estimates for the local density of ISOs currently vary significantly, ranging from approximately 10¹⁴ pc⁻³ (Portegies Zwart et al. 2018) to 8 × 10¹⁴ pc⁻³ (Jewitt et al. 2017), and up to 2 × 10¹⁵ pc⁻³ (Do et al. 2018). Additionally, ’Oumuamua ISSI Team et al. (2019) estimated the underlying mass density of ISOs, considering scenarios involving planetary instabilities. Their results suggest that the total mass of ISOs ejected due to planet–disk interactions ranges from approximately 0.004 to 3 M_⊕,pc⁻³. However, given the current lack of observational data on ISOs, these estimates remain highly uncertain. Nevertheless, the findings presented in these previous studies, as well as those discussed in this work, clearly indicate that planetary systems can eject substantial fractions of their asteroid and comet populations, implying that numerous ISOs populate interstellar space.

We estimate the number of ISOs generated in a single close stellar flyby by combining a size–distribution model of the primordial disk with the ejection fractions f_ej measured in our simulations of single encounters with the compact disk (Sect. 3.2 and Appendix A.1). The latter are defined as the fraction of test particles that become unbound from the Sun, those with e > 1 and positive orbital energy, by the end of each simulation. These values represent the per-encounter efficiency of ejection and scale linearly with the disk mass. We assume a differential size distribution, $$n(D)dD=k{D}^{-q}dD,$$(1)

between diameters D_min and D_max, with bulk density ρ. Here k is a normalization constant that fixes the absolute number of bodies, determined by the total disk mass. For a collisional cascade q ≃ 3.5 (Dohnanyi 1969), close to values inferred in debris disks and the Kuiper belt (Wyatt 2008; Pan & Schlichting 2012). The total disk mass is $$M_{\mathrm{d}}={\int }_{D_{min}}^{D_{max}}\frac{\pi \rho }{6}{D}^{3}n(D)dD=\frac{\pi \rho }{3}k(\sqrt{D_{max}}-\sqrt{D_{min}}),$$(2)

where the prefactor πρ/6 corresponds to the mass of a single spherical body of diameter D. This expression fixes the normalization constant as $$k=\frac{3M_{\mathrm{d}}}{\pi \rho (\sqrt{D_{\mathrm{m}\mathrm{a}\mathrm{x}}}-\sqrt{D_{\mathrm{m}\mathrm{i}\mathrm{n}}})}.$$(3)

With this normalization, the cumulative number of bodies larger than D₀ is $$N(>D_0)={\int }_{D_0}^{D_{max}}n(D)dD\simeq \frac{k}{q-1}{D}_0^{1-q}.$$(4)

This approximation holds for q > 1 and when $D_{\text{max}}^{1-q}\ll D_0^{1-q}$, i.e., when the contribution of the largest objects to the integral is negligible compared to that of D₀ sized bodies.

A stellar flyby that ejects a fraction f_ej = N_ej/N_tot therefore produces $$N_{\mathrm{I}\mathrm{S}\mathrm{O}}(>D_0)=f_{\mathrm{e}\mathrm{j}}N(>D_0).$$(5)

Adopting values of q = 3.5, ρ = 10³ kg m⁻³, D_min = 0.1 km, D_max = 100 km, and D₀ = 1 km, we obtain $$N(>1\mathrm{k}\mathrm{m})\simeq 2.35\times {10}^{11}\left(\frac{M_{\mathrm{d}}}{M_{\oplus }}\right),$$(6)

so that $$N_{\mathrm{I}\mathrm{S}\mathrm{O}}(>1\mathrm{k}\mathrm{m})\simeq 2.35\times {10}^{11}\left(\frac{M_{\mathrm{d}}}{M_{\oplus }}\right)f_{\mathrm{e}\mathrm{j}}.$$(7)

For a 20 M_⊕ primordial disk (e.g., Gomes et al. 2005), we estimate the number of ISOs produced during a single stellar flyby by scaling the ejection fractions f_ej measured from the Compact 1 model. These values correspond to coplanar, prograde encounters at the specified periapsis distances and thus represent upper limits on the ISO yield from individual encounters. $$\begin{array}{rl}& \text{\hspace{0.17em}G-type\hspace{0.17em}}(1M_{\odot },f_{\mathrm{e}\mathrm{j}}=0.024):N_{\text{ISO\hspace{0.17em}}}\sim 1.1\times {10}^{11},v_{\mathrm{\infty }}\simeq 2.5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},\\ & \text{\hspace{0.17em}A-type\hspace{0.17em}}(2M_{\odot },f_{\mathrm{e}\mathrm{j}}=0.025):N_{\text{ISO\hspace{0.17em}}}\sim 1.2\times {10}^{11},v_{\mathrm{\infty }}\simeq 5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},\\ & \text{\hspace{0.17em}B-type\hspace{0.17em}}(5M_{\odot },f_{\mathrm{e}\mathrm{j}}=0.201):N_{\text{ISO\hspace{0.17em}}}\sim 9.5\times {10}^{11},v_{\mathrm{\infty }}\simeq 7\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},\\ & \text{\hspace{0.17em}B-type\hspace{0.17em}}(9M_{\odot },f_{\mathrm{e}\mathrm{j}}=0.409):N_{\text{ISO\hspace{0.17em}}}\sim 1.9\times {10}^{12},v_{\mathrm{\infty }}\simeq 9\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}.\end{array}$$

The resulting asymptotic velocities are comparable to the internal velocity dispersions of young stellar clusters (v_∞ ≃ 2–10 km s⁻¹). This implies that the ejected particles are rapidly mixed into the local stellar environment and contribute efficiently to the interstellar reservoir. Because N_ISO ∝ M_d f_ej, rare but massive perturbers dominate early ISO production by ejecting orders of magnitude more material than solar-type encounters acting on disks of similar mass. Each value of N_ISO in our examples corresponds to the total number of ≳1 km bodies ejected during a single stellar flyby, assuming a 20 M_⊕ primordial solids disk. For instance, a solar-type perturber (1 M_⊙) expels ~10¹¹ objects per encounter, whereas a 9 M_⊙ B-type star can unbind nearly 2 × 10¹² bodies in a single pass. These encounters imply that even a few close interactions during a cluster’s lifetime can substantially enrich the local interstellar medium with cometary debris.

Although in the Compact 1 model the inferred ISO yields per encounter and stellar type are large, they correspond to the total number of kilometer-scale objects ejected during a single, very close stellar flyby (q_⋆ ≲ 300 au). In our simulations, such encounters probe the inner regions or outer edge of the primordial planetesimal disk (see Table 1, Compact 1 model) and occur during the early cluster phase, when disk masses and object numbers are significantly higher than in a fully formed Oort cloud. By contrast, the present-day Oort cloud represents only a small residual fraction of the original planetesimal disk population (Francis 2005). As a result, the number of ISOs produced in an Oort cloud-like structure is expected to be significantly lower. Consequently, the production of ISOs depends sensitively on the encounter impact parameter, the density distribution of particles, and whether the source population is disk-like (early phases) or Oort cloud-like (late phases; e.g., Engelhardt et al. 2017; Pfalzner et al. 2021).

Our results are consistent with theoretical expectations that close encounters and planet–planet scattering efficiently eject particles into interstellar space (e.g., Raymond et al. 2020; Moro-Martín 2018; Torres et al. 2019; Portegies Zwart 2021). Although low-mass stars dominate the IMF, their flybys contribute comparatively little. In contrast, encounters with 5–9 M_⊙ stars, though rare, eject orders of magnitude more bodies, and may thus represent a major source of the first generation of interstellar comets in dense birth clusters.

Fig. 13 Distribution of hyperbolic excess velocities (v∞) of unbound objects from all simulated systems, comparing the Extended (top panel, red) and Compact (bottom panel, blue) models. The histograms represent the particle distributions, while the insets display their cumulative distributions.

5 Summary and conclusions

In this work, we investigated the dynamical evolution of a Solar System analogue (SSA), focusing on the gravitational interactions among four giant planets, a non-self-gravitating particle disk, and passing stars within a dense stellar environment (Sects. 3.1 and 3.2). Our goal was to understand how stellar perturbations reshape the architecture of outer planetary systems, driving the formation of populations analogous to those observed in the outskirts of the Solar System, such as resonant objects, Sednoids, Kuiper Belt and scattered disk objects, Oort cloud-like comets, and ISOs. We explored two disk models: a Compact and an Extended disk, each subjected to either a single close encounter or a sequence of multiple stellar flybys, following the initial conditions summarized in Table 1. The main results of this work are summarized below.

• Disk evolution and cluster environment. In both the Compact and Extended models (Sects. 3.1 and 3.2), the morphology and dynamical heating of the particle disk depend on the encounter geometry and the mass of the perturber. Repeated close encounters with low-mass stars (K and M-type stars) efficiently reproduce structures analogous to the Kuiper belt, scattered disk, and inner Oort cloud (Figs. 3 and 9). These cumulative perturbations drive progressive disk heating and outward migration, raising perihelia and populating long-lived eccentric orbits. Compared to single close encounters, which can strongly reshape or even truncate a disk in one event, the cumulative effect of weaker encounters leads to a gradual redistribution of material rather than its complete removal. As a result, multiple low-mass encounters collectively build up extended scattered populations while preserving the inner regions of the planetary system. For encounters with massive stars, compact disks are comparatively resilient. A-type stars eject only a few per cent of particles, whereas B-type stars truncate the disk at ~20–25 au and excite large eccentricities and inclinations (Fig. A.1). Extended disks, however, are far more fragile; even 2 M_⊙ perturbers remove over half the population, and 5–9 M_⊙ stars can strip 60–75% of bodies while seeding Oort-like orbits (Fig. A.2).

• Encounter geometry and disk populations. We find that the efficiency of disk disruption is strongly geometry-dependent. Coplanar, prograde encounters maximize the transfer of angular momentum and energy to disk particles, whereas inclined or polar encounters are comparatively inefficient (Fig. 7). In the closest encounter we studied (Compact 1 model), a solar-type perturber at q_⋆ ≃ 75 au ejects ~2.4% of particles for i = 0°, giving rise to three distinct dynamical families: a scattered population with low Tisserand parameters (T_U ≲ 2.5) and small perihelia; a scattered tail near the Uranus-scattering boundary (2.5 ≲ T_U ≲ 3.0); and a dominant lower tail that remains dynamically cold at larger perihelia. At higher inclinations (i = 60°−90°), the scattered population nearly disappears and most of the disk survives only mildly perturbed, confirming that the three-branch morphology is a hallmark of prograde or moderately inclined encounters. A similar behavior is found in the Extended 1 models (Sect. 3.1.1). In coplanar encounters, the disk remains largely flattened, whereas polar flybys redistribute angular momentum vertically, producing the isotropic outer populations of an emerging Oort cloud. These results provide a framework for interpreting the architecture of the outer Solar System. When the Sun was still embedded in its birth cluster, repeated stellar encounters with low-mass stars likely drove the outward migration and eccentricity excitation of bodies in its primordial disk, thereby populating the outer regions of the Solar System (see, e.g., Portegies Zwart et al. 2021). Particles scattered onto sufficiently large semimajor axes become susceptible to Galactic tides, which gradually circularize their orbits, while those at smaller a remain on eccentric, long-lived trajectories or become dynamically unstable on secular timescales (e.g., Heisler & Tremaine 1986; Fouchard et al. 2006; Dones et al. 2015; Nesvorný 2018). The tail-like structure observed in semimajor axis–eccentricity space (Figs. 1, 5, 7, and 9) is consistent with this picture, linking the dynamical families formed in our simulations to the seeds of the observed populations of the outer Solar System in its early stages of evolution.

• Implications for the formation of interstellar objects in dense environments. The unbound populations produced in our simulations exhibit characteristic hyperbolic velocities of v_∞ ~ 1–3 km s⁻¹, with a high-velocity tail extending up to ~ 40 km s⁻¹ across both disk models (Fig. 13). These velocities are comparable to the velocity dispersion of stars in open clusters, implying that ejected bodies are rapidly mixed into the local stellar environment. By coupling our kinematic results with the measured ejection fractions from our single-encounter compact-disk simulations (Sect. 3.2.1), we estimated the production of kilometre-sized ISOs during individual stellar flybys (Eq. (5)). For a 20 M_⊕ primordial disk, a solar-type encounter ejects about 10¹⁰ objects per flyby, whereas A- and B-type stars can unbind 10¹¹−10¹² bodies at typical asymptotic velocities of 2–10 km s⁻¹.

  Although these numbers appear large, they represent the total number of bodies ejected from a single planetary system and are therefore strongly diluted when distributed over parsec-scale cluster volumes. The Galactic ISO population must instead arise from the cumulative contribution of planetary systems, mainly from dense stellar environments, where the probability of close encounters is high. Over time, the ISOs will gradually blend into the interstellar medium as clusters dissolve and stars exchange material through mutual encounters (Jílková et al. 2016). Because N_ISO scales directly with disk mass and ejection efficiency, massive stars, though rare, dominate the overall ISO production, while the more frequent flybys of low-mass stars provide a steady but modest background. These results highlight stellar flybys in dense environments as an efficient mechanism for enriching the Galaxy with ISOs.

Our results highlight the key role of stellar encounters in dense environments in shaping the early architecture of planetary systems, seeding distinct dynamical families, and producing ISOs rather than truncating them. This work provides a link between the internal evolution of planetary systems and their external environment. Future studies combining high-resolution N-body simulations with cluster-scale dynamics, stellar evolution, and Galactic dynamics models will be essential for tracing the long-term evolution of SSAs and quantifying how the population of ISOs in dense environments enriches the Galaxy.

Data availability

Movies associated with Figs. 3 and 9 are available at https://www.aanda.org

Acknowledgements

We thank the referee for their suggestions and comments, which helped us improve the quality and clarity of the paper. S.T. thanks Ylva Götberg, Maxwell Cai, Diptajyoti Mukherjee, Simon Portegies Zwart, and Anthony Brown for their valuable feedback and comments. S.T. acknowledges the funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 101034413.

References

Appendix A The effect of A- and B-type stars

In Sect. 3 we examined the dynamical evolution of the SSA under full cluster dynamics. Because of the adopted IMF (see Sect. 2), the vast majority of encounters involved low-mass stars. However, previous studies have shown that massive perturbers can strongly affect the evolution and morphology of planetary systems, and in particular Oort cloud-like structures (e.g., Fouchard et al. 2011). To complement the cluster-based results, we perform dedicated single-encounter simulations with A- and B-type stars, designed to probe the impact of these rare but potentially encounters.

We modeled encounters with stars of 2 M_⊙, 5 M_⊙, and 9 M_⊙ for both the Compact (Appendix A.1) and Extended (Appendix A.2) disk configurations. All simulations adopt coplanar, prograde geometries (i = 0°), which maximize the gravitational interaction between the perturber and the disk (Sect. 3). For consistency with our earlier analysis, we use pericenter distances of 90 au for the Compact 1 model (Sect. 3.2.1) and 300 au for the Extended 1 model (Sect. 3.1.1).

Appendix A.1 The effect of A- and B-type stars in compact disks

In Fig. A.1 we show the outcome of single encounters with A- and B-type stars for the Compact 1 model. This marks the threshold where the 2 M_⊙ A-type star begins to noticeably perturb the disk, ejecting ≃ 2.4% of particles. At larger pericenters, the A-type has little effect, whereas B-type stars at the same distance remove up to ~ 40% of the population. For comparison, at q_⋆ = 100–140 au show negligible excitation, while a 1 M_⊙ perturber at 60 au ejects only ~ 90 particles. Thus, the 90 au case provides a natural reference point for the onset of significant disruption for A- and B-type stars.

As expected, the A-type case (2 M_⊙) is the least disruptive. Around 97.5% of the 2000 particles remain bound, with only 2.5% ejected. The particles in the disk cluster around a_med = 24.1 au (18.4–28.3 au), with small eccentricities (e_med = 0.15) and nearly circular inclinations. About 10% reach e > 0.5, < 0.3% exceed i > 30°, and only a few are ejected beyond the primordial disk (~ 2% to the scattered disk, 0.3% to the inner Oort Cloud). The most distant particles lie at a ≃ 2300 au, while ejected bodies acquire velocities of v_∞ ≃ 5 km s⁻¹.

The 5 M_⊙ encounter is significantly more erosive. 20% of the disk is ejected, with survivors concentrated inside ~ 25 au. The median orbital elements are a_med = 19.9 au (16–84th range 16.4–28.3 au) and e_med = 0.31, with one third of particles reaching e > 0.5. Inclination excitation remains low (i_med ≃ 0°, ~ 1% with i > 30°). About 9.5% of survivors reach the scattered disk and 1.4% the inner Oort Cloud, while the most distant particle extends to ~ 4500 au. The ejected population leaves with v_∞ ≃ 7 km s⁻¹.

The 9 M_⊙ B-type perturber is the most destructive as expected, removing 41% of the disk and producing a sharp truncation at a ≃ 24 au. Bound particles are tightly concentrated around 15–20 au, with a_med = 18.5 au (10.7–33.0 au) and e_med = 0.45. Nearly half of them have e > 0.5, 2.7% exceed i > 30°, and ~ 1% flip to retrograde orbits. The scattered disk and inner Oort Cloud contain 12.5% and 2.2% of the particles, respectively, and the most extreme object reaches a ≃ 1.1 × 10⁵ au. The ejected bodies escape with v_∞ ≃ 9.3 km s⁻¹.

The cumulative ejection fractions F_ej(a₀) shown in the bottom panel of Fig. A.1 confirm that B-type encounters produce systematically stronger stripping at all a₀. The A-type flyby predominantly heats the outer disk (a₀ ≳ 25 au), while the B-type already perturbs material near the inner edge. Despite this, all four giant planets remain bound and nearly circular, underscoring the resilience of the planetary region in encounters with massive stars.

Fig. A.1 Top: Probability density function vs. semimajor axis for bound particles after a single stellar flyby. The histograms are shown for encounters with an A-type star (2 M⊙, blue), a B-type star (5 M⊙, green dashed and 9 M⊙, orange), compared with the initial disk (gray). The inset shows kernel density estimates (KDEs). Bottom: Cumulative ejection fraction as a function of the initial semimajor axis a0 for the same stellar flybys at a pericenter distance of q⋆ = 90 au.

Overall, these results demonstrate the steep dependence of encounter strength on stellar mass. While the A-type flyby leaves the disk largely intact for a very close encounter (90 au), B-type encounters both erode and dynamically excite the population, truncating the disk and seeding long-period orbits. Although only a small fraction of particles reach Oort-like distances in single encounters, such events may nonetheless provide an efficient pathway for populating the early scattered disk and inner Oort Cloud in dense birth environments.

Appendix A.2 The effect of A- and B-type stars in extended disks

In Fig. A.2 we show that encounters with A- and B-type stars are highly disruptive in the Extended 1 model. Between 60% and 76% of the initial particles are ejected into interstellar space. The efficiency of erosion increases systematically with stellar mass, with f_ej = 0.61 for a 2 M_⊙ perturber, f_ej = 0.72 for 5 M_⊙, and f_ej = 0.76 for 9 M_⊙.

Bound particles are confined to a progressively smaller radius with increasing stellar mass. The truncation radius, a₉₀, where 90% of the survivors lie inward, decreases from ~ 250 au (A-type) to ~ 150 au (B5) and ~ 125 au (B9). Median semimajor axes are clustered around a_med = 230–290 au, with long tails extending to tens of thousands of au. Around ~ 1% of the total population, reach a ≃ 2 × 10⁴ au (A-type) and up to ~ 7 × 10⁴ au (B5) or even 3.5 × 10⁴ au (B9). These represent the first seeds of scattered-disk and Oort cloud-like populations.

The eccentricities of the reaming particles in the disk are strongly pumped, with median values e_med ≃ 0.6 in all cases and more than half exceeding e > 0.5. Inclinations are widely dispersed, the A-type flyby produces i_med ≃ 41°, with ~ 58% above 30° and ~ 22% retrograde. The 5 M_⊙ encounter drives i_med ≃ 71° with ~ 66% above 30° and nearly 40% retrograde, while the 9 M_⊙ case pushes the distribution to i_med ≃ 94° and ~ 52% retrograde. For the 2 M_⊙ perturber, 90% of survivors remain inside the initial disk, with only ~ 10% scattered outward (3.2% to 1000 < a < 5000 au, 1.8% to the inner Oort Cloud). In contrast, the 5–9 M_⊙ cases leave only ~ 80% of bound particles in the disk, while ~ 10–15% are placed in the scattered-disk region and up to 5% in the inner Oort Cloud. No significant population reaches the outer Oort Cloud in these single encounters.

The ejected population leaves with characteristic asymptotic speeds of v_∞ ≃ 2–3 km s⁻¹, somewhat lower than in the compact-disk case (Appendix A.1). This difference reflects the weaker binding of the extended disk: particles at large radii require only modest perturbations to escape, yielding high ejection fractions but relatively low excess velocities. These results show that single close encounters with massive stars (A- and B-type) are a highly efficient channel for producing interstellar comets from extended disks at very close passages (~ 300 au). While A-type stars already eject more than half the population, B-type encounters both erode the disk more strongly and drive more extreme excitation, seeding scattered-disk and Oort-like reservoirs and yielding 10¹¹−10¹² ISO-scale bodies for a 20 M_⊕ disk (Sect. 4).

Fig. A.2 Top: Probability density function vs. semimajor axis for bound particles after a single stellar flyby for the Extended 1 model. The histograms are shown for encounters with an A-type star (2 M⊙, blue), a B-type star (5 M⊙, green dashed and 9 M⊙, orange), compared with the initial disk (gray). The inset shows kernel density estimates (KDEs). Bottom: Cumulative ejection fraction as a function of the initial semimajor axis a0 for the same stellar flybys at a pericenter distance of q⋆ = 300 au.

All Tables

All Figures

	Fig. 1 Extended 1 model. Semimajor axis as a function of eccentricity and particle count after an encounter with a 1 M⊙ star at 300 au. The color bar indicates the perihelion distance of the particles. Each panel corresponds to a different encounter inclination angle (0°, 30°, 60°, and 90°).
In the text

Fig. 2 Extended 1 model. Cumulative distribution of the final energy of particles in the disk. The colored lines correspond to different inclination angles of the encounter (0°, 30°, 60°, and 90°), while the black dotted line represents the initial particle distribution as defined in Table 1. Negative energy values indicate particles that remain bound to the system after the encounter, whereas positive values correspond to ISOs. The Kolmogorov-Smirnov (KS) probabilities relative to the initial distribution for the different angles are 0.90% for 0°, 0.19% for 30°, 0.07% for 60°, and 0.05% for 90°.

In the text

Fig. 3 Extended N model. Semimajor axis as a function of eccentricity (top panel) and orbital inclination (bottom panel) for all the particles in the 200 simulated systems. The gray areas represent the different regions of the Solar System: Kuiper belt (KB), scattered disk (SD), Sednoids, inner and outer Oort cloud (OC); the dashed lines show the 2:3, 1:2, and 2:5 resonances with Neptune. The integration time is set to 100 Myr. An animation is available online.

In the text

	Fig. 4 Extended N model. Orbital evolution of particles in the disk over 100 Myr. The bottom panels in each plot depict the semimajor axis as a function of time, while the top panels show the distance of the perturber. The panels represent systems 128, 14, 100, 86, and 157 (as labeled). The colored lines correspond to individual particles within each system.
In the text

Fig. 5 Extended N model. Orbital elements of particles in the disk after 100 Myr and multiple stellar encounters. First column: semimajor axis as a function of eccentricity, with particles color-coded by their perihelion distances. The green, blue, and gray shaded regions highlight different populations formed due to stellar encounters. Second column: perihelion as a function of orbital inclination, with dots color-coded by aphelion distance. The gray shaded areas represent the distinct regions of the Solar System: Kuiper Belt (KB), scattered disk (SD), Sednoids, and the inner Oort Cloud (OC). Third column: distribution of orbital energy for the particles. The red histograms correspond to the initial energy distribution, and the blue curve represents the final energy distribution. The blue curves with positive values indicate interstellar comets. Each row corresponds to systems numbered 128, 14, 100, 86, and 157, respectively.

In the text

	Fig. 6 Extended N model. Phase-space distribution (X vs. Y) of particles in the disk after 100 Myr. The color scale indicates the particles’ positions along the vertical (Z) axis. The panels (from top to bottom) show systems 128, 14, 100, 86, and 157.
In the text

	Fig. 7 Compact 1 model. Semimajor axis as a function of eccentricity after an encounter with a 1 M⊙ star at q⋆ ≃ 75 au. The color bar indicates the perihelion distance of the particles. Each panel corresponds to a different inclination angle of the encounter (0°, 30°, 60°, and 90°).
In the text

	Fig. 8 Compact 1 model. Cumulative distribution of the final orbital energy of disk particles. The colored lines correspond to different inclination angles of the encounter (0°, 30°, 60°, and 90°), while the black dashed line represents the initial distribution. Negative values (−1/a < 0) indicate bound orbits, while positive values correspond to ejected ISOs.
In the text

Fig. 9 Compact N model. Semimajor axis as a function of eccentricity (top panel) and orbital inclination (bottom panel) for all the particles in the 200 simulated systems. The gray areas represent the different regions of the Solar System: Kuiper belt (KB), scattered disk (SD), and the Sednoids. The dashed lines show the mean motion resonances with Neptune. The integration time is set to 100 Myr. An animation is available online.

In the text

	Fig. 10 Compact N model. Orbital evolution of particles in the disk over 100 Myr. The bottom panels in each plot depict the semimajor axis as a function of time, while the top panels show the distance of the perturber. The panels represent systems 128, 14, 100, 86, and 157 (as labeled). The colored lines correspond to the individual particles within each system.
In the text

Fig. 11 Compact N model. Orbital elements of particles in the disk after 100 Myr and multiple stellar encounters. First column: semimajor axis as a function of eccentricity. The particles are color-coded by their perihelion distances. The green, blue, and gray shaded regions highlight different populations formed due to stellar encounters. Second column: perihelion as a function of orbital inclination. The dots are color-coded by aphelion distance. The gray shaded areas represent distinct regions of the Solar System: Kuiper Belt (KB) and scattered disk (SD). Third column: distribution of orbital energy for the particles. The red histograms correspond to the initial energy distribution, and the blue curve represents the final energy distribution. The blue curves with positive values indicate interstellar comets. The rows (from top to bottom) correspond to systems 128, 14, 100, 86, and 157.

In the text

	Fig. 12 Compact N model. Phase-space distribution (X vs. Y) of particles in the disk after 100 Myr. The color scale indicates the particles’ positions along the vertical (Z) axis. The panels (from top to bottom) show systems 128, 14, 100, 86, and 157.
In the text

	Fig. 13 Distribution of hyperbolic excess velocities (v∞) of unbound objects from all simulated systems, comparing the Extended (top panel, red) and Compact (bottom panel, blue) models. The histograms represent the particle distributions, while the insets display their cumulative distributions.
In the text

Fig. A.1 Top: Probability density function vs. semimajor axis for bound particles after a single stellar flyby. The histograms are shown for encounters with an A-type star (2 M⊙, blue), a B-type star (5 M⊙, green dashed and 9 M⊙, orange), compared with the initial disk (gray). The inset shows kernel density estimates (KDEs). Bottom: Cumulative ejection fraction as a function of the initial semimajor axis a0 for the same stellar flybys at a pericenter distance of q⋆ = 90 au.

In the text

Fig. A.2 Top: Probability density function vs. semimajor axis for bound particles after a single stellar flyby for the Extended 1 model. The histograms are shown for encounters with an A-type star (2 M⊙, blue), a B-type star (5 M⊙, green dashed and 9 M⊙, orange), compared with the initial disk (gray). The inset shows kernel density estimates (KDEs). Bottom: Cumulative ejection fraction as a function of the initial semimajor axis a0 for the same stellar flybys at a pericenter distance of q⋆ = 300 au.

In the text