© The Authors 2025

1 Introduction

Planetary systems form within circumstellar disks mainly composed of gas and a minor fraction of dust. Gravitational interactions between the gaseous disk and embedded planets drive planetary migration. The direction and speed of this migration depend on the planet’s mass and the local properties of the gaseous disk (Goldreich & Tremaine 1979; Lin & Papaloizou 1986).

When the relative migration rate is convergent and adiabatic (with the orbital period ratio decreasing over a timescale longer than the resonant libration period), planets in low-eccentricity orbits are expected to become trapped in mean-motion resonances (MMRs) (Baruteau et al. 2014; Nelson 2018). When more than two planets are formed in the same system, a resonant chain is expected to form, capturing all adjacent pairs in resonance¹.

However, out of ~5500 known planetary systems, around 1000 host multiple planets, yet fewer than 30 exhibit resonant chains. Figure 1 illustrates systems near first-order resonances (only those within 5% of period ratios implying 2:1 and 3:2 MMRs), based on data retrieved from the Exoplanet Catalog (https://exoplanet.eu/ as of November 2024).

The data reveal that, with the exception of HR 8799, all systems exhibiting resonant chains consist of super-Earths with orbital periods typically within 200 days of their host star. In contrast, cold Jupiter systems – a population of giant planets with semi-major axes beyond 1 au – are frequently observed in isolation (or at least, lacking close companions, e.g., with period ratios below ~5, Zhu 2022; Rosenthal et al. 2024). Their significant eccentricities (median of 0.2, e.g., Kipping 2013; Bowler et al. 2020; Rosenthal et al. 2021) largely suggests that planetary ejections from planet-planet scattering has shaped their orbits and decreased their multiplicities (see Weidenschilling & Marzari 1996; Chatterjee et al. 2008; Juric’ & Tremaine 2008; Davies et al. 2014, and references therein).

Yet, a small fraction (15%) of these cold Jupiters are observed alongside other cold Jupiter companions, frequently observed near the 2:1 MMR and to a lesser extent near the 3:2 (Charalambous et al. 2022; Cerioni et al. 2022). This leads to the following central questions: what mechanisms can disrupt long-period resonant chains of systems with giant planets? Could isolated cold Jupiters result from the breakup of such chains?

For close-in super-Earth planets, several processes have been identified that prevent the formation of resonant pairs, such as overstable librations (Goldreich & Schlichting 2014), disk turbulence (Batygin & Adams 2017), planetesimal-driven migration (Chatterjee et al. 2008), and tidally driven resonance divergence induced by the disk (Choksi & Chiang 2020) or tides acting on the planets (Lithwick & Wu 2012; Batygin & Morbidelli 2013; Charalambous et al. 2023). Recently, Dai et al. (2024) showed that resonant configurations are practically universal among young planetary systems. Over time, however, these resonances tend to dissolve. This supports the idea that planets are initially assembled in resonant configurations, which are gradually disrupted by various processes, including dynamical instabilities such as those suggested by the “breaking-the-chains” model (Izidoro et al. 2017; Goldberg & Batygin 2022).

For cold Jupiters, the mechanisms responsible for breaking MMRs remain more elusive. Tides are inefficient (Lithwick & Wu 2012; Batygin & Morbidelli 2013), and planetesimal-driven migration and turbulence are less effective due to the higher planetary masses (Chatterjee et al. 2008). One possibility is that planets never got captured in the first place. However, the recent numerical experiments by Nagpal et al. (2024) show that the conditions to disrupt the resonances for giant planets are somewhat stringent, demanding fast migration rates and inefficient eccentricity damping rates² – conditions not naturally accommodated within the type-II migration paradigm. Furthermore, in systems that still retain a disk, such as PDS 70 and HD 163296, observational evidence suggests that long-period resonant chains can form (Bae et al. 2019; Garrido-Deutelmoser et al. 2023), while young systems such as HR 8799 still preserve such resonance structure (Goz’dziewski & Migaszewski 2020; Zurlo et al. 2022). As such, resonance chains for cold Jupiters may indeed be broken after the disk dispersal phase by a process that has yet to be determined.

In this work, we propose a mechanism analogous to the breaking-the-chains model for gas giant planets, where resonant chains are formed through disk-driven migration. Our key contribution is to consider the system as non-isolated, evaluating the role of stellar flybys in breaking these resonances.

Most models assume that stars are isolated and that the process of planet formation and evolution occurs almost unperturbed. However, it is well known that following the collapse of molecular clouds, stars often form multiple stellar systems (Pringle 1989; Reipurth et al. 2014; Offner et al. 2023), which evolve within stellar clusters and associations (Lada & Lada 2003). In these regions, stellar encounters are not only possible but rather frequent. The likelihood of these events strongly depends on the adopted flyby definition and it increases with higher stellar density (Pfalzner 2013; Winter et al. 2018), as expected. Here, we follow the typical flyby definition formulated in Cuello et al. (2023): a stellar encounter occurs whenever two stars – with eccentricity equal to or greater than 1 – get closer than 1000 au. For instance, for a 3 Myr old solar-like star in a leaky cluster or in an OB association, the probability of experiencing a flyby is around 20% according to Pfalzner (2013). By that time, the process of planet formation should already have taken place (Manara et al. 2023), meaning that stellar flybys can impact fully formed planetary systems, which may already be in resonance.

Scenarios where a mature planetary system is perturbed by an intruding star have been studied over the years. Laughlin & Adams (1998) showed that stellar encounters in dense open clusters can cause significant orbital disruption, which could explain the presence of Jovian planets in eccentric and relatively close orbits. In a similar vein, Malmberg et al. (2011) explored through numerical simulations how stellar flybys can trigger scattering events among planets, disrupting stable configurations. Li & Adams (2015) compare the effects of single and binary stars, finding that binaries generally have a greater impact on planetary systems, particularly those with widely spaced planets. Li et al. (2019, 2020) further explore scenarios in which planetary systems lose planets due to catastrophic flybys, where the passing star, that has its own planetary system, is allowed to cross between the planets of the affected system, significantly disrupting their orbits.

Our approach is novel since it focuses on mature systems with resonant chains in wide orbits, while previous studies examined compact systems in their birth clusters. Here, we explore systems in the field with giant planets perturbed by a single passing star, which, despite fast relative velocities, has a non-negligible impact (Veras & Moeckel 2012). Additionally, we analyze the effects of the passing star at realistic distances in different resonant chains and assess their long-term fates. The remainder of this paper is organized as follows. Section 2 outlines the typical flyby parameters, and details the setup for the simulations. In Section 3, we discuss the results of our simulations after varying several parameters, and in Section 4 we summarize our findings.

Fig. 1 Systems within 5% of resonant chains involving the 3:2 and 2:1 period ratio commensurabilities between adjacent planet pairs. Data have been taken from the Exoplanet Catalog (https://exoplanet.eu/catalog/). The size of the dots indicates planetary masses: large dots represent planets with masses of at least 1 mJup, while smaller dots correspond to less massive planets.

2 Flyby parameters and setup of the simulations

We consider systems of three planets on coplanar orbits around a central star with mass M_* = 1 M_⊙. We denote the inner, middle, and outer planetary masses as m₁, m₂, and m₃, respectively, with M_* ≫ m₁, m₂, m₃. The three planets in the system are inter-connected through resonant chains: planets 1 and 2 are in a mean-motion resonance, as are planets 2 and 3. Consequently, the entire system exhibits a three-planet resonance, a multi-resonant state. Our focus is primarily on zero-order three-planet resonances (q_3pl = 0, also known as Laplace resonance), as these are the only ones observed among exoplanets at present³.

Since we are interested on the effects of the flyby on the system’s dynamics, we force the planets into a resonant chain, mimicking in situ formation and only experiencing short-term migration. We adiabatically push the outer planet inward allowing it to capture the inner planets in resonance, leading to the libration of the two- and three-planet resonant angles. Orbital and eccentricity damping timescales of τ_a = 10⁶ yrs and τ_a/τ_e = 100 are applied to account for planet-disk interactions, with the planets initially placed 10% away from resonance. The capture takes place at the beginning of the simulation, after which we let the system evolve in isolation for 5 × 10⁵ yrs, influenced only by their gravitational interactions to ensure stability. After this period, a passing star with a mass of M_fb = 1 M approaches the system. It is important to mention that in the absence of a stellar flyby, both the 2:1-2:1 and 3:2-3:2 configurations remain stable and locked in resonance for at least 10⁷ orbits of the innermost planet, ~5 × 10⁸ years, with typical equilibrium eccentricities lower than 5%, and the libration amplitude of the resonant angles around 20°.

We consider that a flyby occurs whenever an unbound star from the field gets closer than 1000 au to the planet-hosting star (Cuello et al. 2023). In order to avoid spurious gravitational effects on the planetary system under study and following Bhandare et al. (2016), we initialize the perturber star at 20 times the distance of minimum approach (also referred to as orbital pericenter). Breslau et al. (2017) showed that considering different azimuthal positions for a given planet when the perturber goes through pericenter leads to different dynamical outcomes (perturbation, ejection or capture). To include this effect, for a given numerical setup, we perform eight 45-degree rotations of the initial true anomaly (position along the orbit). This approach allows us to obtain meaningful results since, for a given resonant layout, we take into account how different initial conditions affect the dynamical outcome.

To sample a limited but representative region of the parameter space, we consider individual stellar perturbers with masses and pericenter distances ranging from 0.5 to 2 M_⊕ and from 75 to 500 au, respectively. These values are consistent with recent calculations of stellar cluster dynamics (Schoettler & Owen 2024). For instance, Schoettler et al. (2020) find that more than half of the encounters occur at distances closer than 300 au in the Orion Nebula Cluster. In more sparse environments, stellar encounters still occur at a moderate rate with low velocities and at pericenter distances of several hundred au (Pfalzner 2013; Winter et al. 2018, 2024). It is worth noting that the likelihood of an encounter decreases with the cluster age. However, as will be shown in this work, under certain conditions, a distant flyby is enough to disrupt a resonant chain of planets. If the system would contain gas, the perturbative effects would be damped (e.g., Marzari & Picogna 2013; Nealon et al. 2020). This is why we decide to focus on gas-less planetary systems. Therefore, the effects reported here should be considered as upper limits for planetary systems embedded in gaseous disks. In addition, as done in Cuello et al. (2019), we assume that there is no preferential direction in which the flyby occurs, and we therefore consider 5 possible orbital inclinations with respect to the initial planetary system orbital mid plane: (i) coplanar prograde, (ii) inclined prograde, (iii) polar or orthogonal, (iv) inclined retrograde, and (v) coplanar retrograde.

In this scenario, the velocity at which the perturber crosses the sphere of 1000 au centered around the planet-hosting star determines the type of encounter. For very high speeds, one can use the impulse approximation assuming that the planets do not progress in their orbits around the central star during the encounter. Typically, these flybys are penetrating or catastrophic. On the contrary, if the perturber has a large pericenter distance and arrives at a lower speed, then we can describe the perturbation in the secular regime as the planets complete various orbits during the interaction. In such a case, the flyby is expected to lead to variations in the orbital eccentricity while the planetary semi-major axis should remain almost unchanged (unless the perturber crosses the planetary orbits on a parabolic orbit).

Taking into account all the points mentioned above, our fiducial simulation involves a system consisting in three Jovian planets with (0.75, 1, 1.25) m_Jup in a resonant chain, and a one solar mass star approaching the planetary system with different inclinations β ∈ [0°, 45°, 90°, 135°, 180°]. For the minimum distance between the stars, we consider six cases: q ∈ [75, 100, 125, 150, 175, 200] au. We analyze two specific first-order resonance chains, namely the 3:2-3:2, and 2:1-2:1. In each setup, the outer planet is positioned at ~40 au, ensuring that the geometric cross section of both layouts remains the same. With these parameters, the encounter can be regarded as a secular perturbation (see Heggie & Rasio 1996; Spurzem et al. 2009). For each combination, we perform eight rotations per system modifying the angular position of the planets along their respective orbits, resulting in a total of 480 simulations. We perform our simulations using the N-body code REBOUND (Rein & Liu 2012) employing the IAS15 integrator (Rein & Spiegel 2015), incorporating a stellar flyby on a parabolic orbit. Once the flyby is over (i.e., when the perturber is at more than 1000 au from the primary star), we remove the perturber from the simulation. Finally, we evolve the planetary system for another 10⁵ orbits of the innermost planet, using the hybrid integrator TRACE (Lu et al. 2024). TRACE uses the symplectic WHFast integrator (Rein & Tamayo 2015) when particles are well separated, and automatically switches over to a Bulirsch-Stoer/IAS15 integrator during close encounters.

3 Results

Given that giant planets on wide orbits around solar type stars may not be uncommon (Fernandes et al. 2019; Fulton et al. 2021), our working hypothesis is that these planets initially formed in resonant chains on wide orbits. In this section, we examine the possibility that the system maintains its resonant configuration after an interaction with a nearby star. Section 3.1 focuses on the evolution within the 2:1-2:1 and 3:2-3:2 resonant chains, while Sections 3.2 and 3.3 discuss the impact of planetary and stellar masses on these interactions.

3.1 Effects of a stellar flyby on wide orbit resonant chains

3.1.1 The immediate response

In the top panels of Figure 2, we present the outcomes of our simulations for systems initially trapped in a 3:2-3:2 (left) and a 2:1-2:1 resonant chain (right). The color coding represents the survival probability of period ratios, with darker colors indicating higher probabilities and lighter colors indicating lower probabilities. Crosses denote instances of ejections or collisions with the star, resulting from destabilization due to the flyby. The y-axis in both panels shows the inclination of the passing star, while the x-axis shows the distance of minimum approach normalized by the position of the last planet in the system (40 au in all cases).

We assess the systems’ stability by considering whether they maintain their orbital periods within 5% of the fixed point after the passing star moves beyond a distance of 1000 au from the central star, a criterion we refer to as the immediate response. Although the outer planet in both scenarios is placed at the same distance and, therefore, similarly perturbed by the passing star, Figure 2 shows that the 3:2-3:2 layout is more susceptible to disruption than the 2:1-2:1. In the 3:2-3:2 setup, the planets are closer to each other when the flyby occurs. This closer proximity, combined with the compact nature of these configuration, makes the system more vulnerable to destabilization by stochastic perturbations.

When the passing star approaches the system in the same direction as the planets’ orbital motion (β = 0°), it disrupts the resonant chain more effectively than during a retrograde approach, even at greater encounter distances. Li et al. (2019) observed a similar effect, where prograde encounters were more likely to destabilize the system and allow the intruding star to capture planets. By integrating across all inclination values, we obtain the cumulative distribution distribution (CDF) of q/a₃ for all possible outcomes: collision and ejection, preservation of orbital periods, and overall stability, as shown in the bottom panels. The immediate response of the system to the effects of the flyby in the 2:1-2:1 resonant chain shows that whenever q/a₃ exceeds 3.5, the systems maintain their constant period ratios, while in the 3:2-3:2 configuration, the systems retain the resonant relation for q/a₃ > 4.5.

As mentioned in Section 1, the fact that the period ratios of consecutive pairs lay in (or close to) an integer ratio does not necessarily mean the system is resonant. For that purpose, the resonant angles must librate around a fixed value. In Figure 3 we show a typical temporal evolution of the planets immediately after the close encounter with the passing star. In this example, the pericentric distance of the approaching star is q = 125 au (q/a₃ ~ 3.5), with an inclination of β = 45°. As expected, the outermost planet is most affected by the interaction with the stellar flyby, causing a bigger instantaneous change in P₃/P₂ than in P₂/P₁. However, while the period ratios remain locked in a 2:1-2:1 commensurability after the passage of the star, the resonant angles begin to circulate, meaning that the system has, in fact, left the resonance.

In such cases, the systems are left with commensurate period ratios but lack the protection provided by the resonance (see Chatterjee et al. 2008). The final orbital elements after the flyby differ significantly from those just after the resonant capture, which jeopardizes the system’s dynamical stability. While secular timescales are orders of magnitude longer than orbital timescales, in this first set of simulations, we do not integrate over such extended periods. Instead, we keep track of the resonant angles for several libration cycles following the encounter. However, to assess the system’s ability to maintain constant period ratios and evaluate their long-term outcomes, it is essential to consider the secular evolution. Therefore, we conduct additional simulations, extending the integration by 10⁵ orbits of the innermost planet.

Fig. 2 Top: survival probability of period ratios (shading from lighter to darker) immediately after a star with a mass of Mfb = 1 M⊙ flies by a three-planet system with (0.75, 1, 1.25) Jupiter masses in 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains and leaves the surroundings of the system to a distance >1000 au, considering different inclinations and distances relative to the orbital plane of the planets. Bottom: cumulative distribution function (CDF) for the possible outcomes of the simulations: collisions with the star and ejections (red dashed line), survival of the period ratios (their resonant angles not necessarily librating, in solid purple), systems where period ratios do not remain commensurate (green), plotted against the pericentric distance of the flyby star relative to the outermost planet.

Fig. 3 Typical time evolution for the three planets in the 2:1-2:1 resonant chain after a close encounter with a flying-by star with mass 1 M⊙ and inclination β = 45°. The distance of minimum approach between the two stars is q = 125 au. In the top panels from left to right, we show the period ratio between the adjacent pairs, the eccentricities and inclinations, while in the bottom panels, the two-body resonant angles for each pair (left and middle plots) and the three-body resonant angle (in the bottom right panel).

Fig. 4 Delayed eccentricities (top row), inclinations (middle), and semi-major axis (bottom row) for all simulations in the 3:2-3:2 (left column) and 2:1-2:1 (right column) resonant chains. Colors indicate the pericentric distance q between the stars during the flyby event, with each color corresponding to a specific q, while the symbols represent the inclination β of the flyby star. The black dashed line indicates the stability limit given by the orbit-crossing value ecross.

3.1.2 The delayed response

Figures 4 and 5 present the results of the delayed interactions between the passing star and the planets, with pericentric distance of the flyby represented in different colors. The symbols show the flyby’s inclination, and the transparency of the symbols depicts the planets’ positions in the system prior to the encounter. From top to bottom of Fig. 4, we show the eccentricities, inclinations and semi-major axis, respectively. For clarity, we exclude planets that end the simulations beyond 400 au, which is 10 times the initial distance of the outermost planet. In Fig. 5 we show the systems that survive the flyby in the period ratio plane.

For the chosen parameters, the encounter between the passing star and the planetary system is slow and tidal (q ≫ a₃), occurring within the secular regime. For a single planet, the change in the eccentricity can then be expressed using the formulation by Heggie & Rasio (1996) as $$\delta e=-\frac{15\pi }{16}{\left(\frac{2M_{\text{fb}}^2{a}^3}{M_{123}{M}_{12}{q}^3}\right)}^{1/2}e\sqrt{1-e^2}\sin (2\mathrm{\Omega }){\sin }^2\beta ,$$(1) where M₁₂ = M_* + m_pl and M₁₂₃ = M_fb + M₁₂. Encounters at larger relative distances cause smaller changes in eccentricity, while closer encounters result in more significant perturbations. This power-law dependence on the ratio (a/q)^3/2, which we have verified for a single planet through simulations, shows that even distant stellar flybys, despite being slow, can induce long-term changes in the orbital parameters of planetary systems. The strength of the perturbation decreases with increasing distance from the passing star, making closer encounters far more disruptive to the planet’s orbit.

In the case of systems with three planets, these assumptions still hold: the encounter is slow and tidal relative to the outermost planet, and the perturbations happens in the secular limit. As shown in Figure 4, the changes in semi-major axis are less significant than those in eccentricities, which is the most perturbed orbital element. The important thing here is the amount of energy exchanged during the encounter, even very small changes in eccentricity can have profound effects, potentially leading to close encounters and, eventually, the disruption of the system. The closer the encounter, the more significant the disruptions to the orbits of any planets in the system. Such events may have played a major role in shaping the currently observed exo-planet population, particularly by contributing to the scarcity of wide-orbit systems in 3:2-3:2 resonant chains.

In all simulations where eccentricities exceed the crossing-orbit value, e_i > e_cross ≃ 2k/3 j for a j : j - k MMR (e.g., Tamayo et al. 2021), planets are eventually ejected, ultimately breaking the chains or leaving one planet at a distance dynamically decoupled from the inner pair. In the 3:2-3:2 chain, these changes are particularly pronounced, even during distant flybys, making it more difficult to find exoplanetary systems in this configuration on wide orbits, even if the event is not catastrophic.

Scattering within the systems triggers subsequent orbital instabilities, where planets gain eccentricity, leading to crossing orbits (simulations that end above the dashed lines in the top panels). This process often results in planets either being ejected or colliding with the host star, even long after the stellar perturber has left the planetary system surroundings. Closer encounters, as well as more distant encounters in prograde orbits, induce more significant changes in the orbital elements.

It is important to note that while Figure 4 does not show simulations with flyby pericentric distances beyond q > 200 au, we conducted simulations with flybys up to q = 1000 au. Fly-bys at such large distances should exert negligible gravitational effects on planetary systems, with typical sizes between 10 and 100 au. However, in the 3:2-3:2 resonance setup, we find that although systems with q/a₃ > 4.5 (q = 175 au) initially survive the encounter, only some remain stable over secular timescales. This results from the marginal stability of such massive planets in compact configurations, where even minor perturbations can break the system. We also observe that, although not shown here, reducing the planetary masses improves the stability of this configuration. A Jupiter-mass system in the 3:2-3:2 resonant chain is at the brink of instability, and the flyby provides just enough perturbation in the eccentricities (less than 1%) to disrupt the stability provided by the resonance.

In contrast, nearly all period ratios in the 2:1-2:1 chain remain within 5% of the exact integer value once the pericentric distance of the passing star exceeds q > 125 au (equivalently q/a₃ ≳ 3.5). However, as with the 3:2-3:2 chain, the resonant angles leave the commensurability. Although the value we obtain for q/a₃ is slightly smaller than that reported in Li & Adams (2015), who examined four-planet systems perturbed by a passing binary and a direct comparison is not possible, their observation that the cross-sections for planetary systems interacting with single stars are smaller than those for binary encounters suggests that our results are consistent with theirs. Besides, note that for q = 125 au, when the flyby is coplanar (represented by circles), the planets are left on inclined orbits. This contrasts with the case of a star passing near a protoplanetary disk, where, after relaxation, the particles return to the plane of the disk (Marzari & Picogna 2013; Picogna & Marzari 2014; Nealon et al. 2020).

In the cases in which the flyby disrupts the chain and the systems go unstable, the effect for the two different configurations generates different eccentricities for the remaining planets. The eccentricities acquired by planets that survive the extended simulations are higher if the system formed near the 3:2-3:2 resonant chain compared to those originally trapped in the 2:1-2:1.

When comparing the immediate versus delayed effects of the flyby (Figures 2 and 4, respectively), we find little difference in the 2:1-2:1 resonant chain, but observe a significant difference in the 3:2-3:2 setup. Analyzing systems within the 3:2-3:2 resonance immediately after the star departs might suggest that the planets retain their period ratios. However, over time, planet-planet interactions eventually destabilize the system. As previously mentioned, although the encounter between the passing star and the planetary system is secular, the immediate effect of the flyby is most noticeable in eccentricities, which induces long-term variations in the semi-major axes. In contrast, the long-term effects on the 2:1-2:1 resonance seem to have minimal impact on the overall configuration, as shown in Figure 5. However, while the period ratio and orbital elements remain largely unaffected, the resonant angles circulate.

In the work of Cerioni et al. (2022), it is noted that large-mass systems show a strong preference for orbital period ratios close to the 2:1 commensurabilities, rather than for the 3:2 and more compact configurations. Since the long-term effects of the flyby are considerably stronger in the 3:2-3:2 resonant chain compared to the 2:1-2:1, particularly in producing substantial changes in eccentricities, it suggests that resonant chains in the 3:2-3:2 layout are unlikely to survive a flyby. Therefore, systems in this configuration are less likely to be observed, especially considering the high probability of a past close encounter between stars (Pfalzner 2013).

To describe the dynamics of a single zeroth-order three-planet resonance, it is possible to write an integrable one-degree-of-freedom Hamiltonian to second order in eccentricities (Quillen 2011; Charalambous et al. 2018; Petit et al. 2020). This formulation uses the three-planet angle ϕ_3pl (see Section 1), which is independent of ϖ_i. As a result, the system preserves its angular momentum deficit (Laskar 1997). Consequently, if the system is only influenced by these resonances, initially circular orbits will remain circular.

When the resonances are well isolated, there is no possibility for large-scale chaos. In such cases, the system behaves as nearly secular and can be considered as long-term stable. However, due to the perturbation introduced by the passing star, the eccentricities get excited, and the systems can no longer be considered stable (see Figure 4). In the absence of a stellar flyby, the systems in our simulations remain stable for at least 10⁸ years, with all angles librating. However, the flyby causes the three-planet resonant angle to circulate, triggering changes in eccentricities and leading to the observed instabilities in the 3:2-3:2 resonant chain. Opposite to this, in the 2:1-2:1 configuration, even with the Laplace angle circulating, the systems are able to maintain the eccentricities low enough that instabilities do not arise.

In summary, while the immediate response of planetary systems to a nearby stellar flyby may appear similar across the configurations analyzed, the long-term outcomes differ significantly. Systems initially near a 2:1-2:1 resonant chain are more likely to retain their period ratios, implying a resonant configuration after long-term interactions, albeit often with circulating resonant angles. In contrast, systems starting near a 3:2-3:2 configuration experience planet-planet scattering triggered by the flyby, producing a population of eccentric cold Jupiters from orbits that would otherwise remain stable. Lastly, recent experiments by Maas et al., (in prep.) including a sequence of flybys from cluster simulations applied to HR 8799 showed that in the absence of resonant chains (i.e., with non-librating resonant angles), the survival probability of HR 8799-like systems is nearly zero, a stark contrast to the higher stability observed in resonant systems.

Fig. 5 Conditions leading to configurations within 5% of the 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains in the period-ratio plane (P2/P1, P3/P2). Relevant two- and three-planet MMRs are indicated by dashed gray curves. Same colors and markers as in Fig. 4.

Fig. 6 Top panels: Period ratio survival for sub-Saturn mass planets, in 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains. Bottom panels: delayed survival (after 105 orbits of the innermost planet) in the period-ratio plane. Colors and markers follow the scheme used in Figure 4.

3.2 The effect of the mass of the planets

3.2.1 Sub-Saturn mass planets

In Section 3.1, we examined the dynamical evolution of the system post stellar flyby, focusing on scenarios with Jupiter-mass planets in different resonant configurations. In this section, we extend our analysis by investigating how varying the individual planet masses impacts the stability of the system. To achieve this, we replicate our fiducial set of simulations, but instead of Jupiter-mass planets, we consider sub-Saturn mass planets of (40, 53.33, 71.11) m_⊕. While the mass ratios between the planets are kept constant, the ratios between the planets and the central star have changed.

In the top panels of Figure 6, we present the immediate response of the planetary orbits to the stellar flyby. On the left, we show the case of a 3:2-3:2 resonant chain, while the right displays the 2:1-2:1 resonance. The bottom panels depict the delayed response of the system in the period ratio plane. Consistent with our fiducial Jupiter-mass simulations, a retrograde stellar flyby continues to enhance the likelihood of the planetary system’s to maintain constant period ratios. Opposite to what we observe for giant planets, for systems of lower-mass planets, flybys beyond q = 125 au still permit the system to remain intact with period ratios within 5% of the 3:2-3:2 configuration, at least in the considered timescales, although the resonant angles are not librating. This is in contrast to the Jupiter-mass case, where more distant flybys destabilize the system. Collisions the star and ejections happen when the star passes closer to the system compared to the case of the giants.

Additionally, for closer encounters, the final configuration of the planets shows a broader dispersion in the period ratios P₂/P₁ and in P₃/P₂, centered around the location of the three-planet MMR, particularly in the 2:1-2:1 case. This spread occurs independently of the inclination of the flyby and highlights the greater robustness of systems with smaller planets to close stellar encounters.

In short, lower-mass planets exhibit a higher degree of resilience to perturbations induced by stellar flybys compared to their giant counterparts. This suggests that a multi-resonant state involving low-mass planets can survive in systems where stellar encounters occur at relatively large distances, opening the possibility of finding such systems in wide orbits.

3.2.2 Super-Jupiter mass planets

Finally, we perform a last set of simulations considering higher mass systems, keeping fixed mass ratios between the planets. For the rest of the Section, we continue the analysis and comparisons for the 2:1-2:1 resonant chain, since we are mostly interested in giant planet systems that may survive the stellar flyby, which is not the case of the 3:2-3:2 configuration in this mass range.

We explore scenarios where the fiducial planetary masses are scaled by factors of 2, 3, 4, and 5, resulting in systems with (1.5, 2, 2.5) m_Jup, (2.25, 3, 3.75) m_Jup, (3, 4, 5) m_Jup and (3.75, 5, 6.25) m_Jup. The effects of the delayed response of these higher-mass planets to the stellar flyby are shown in Figure 7, where we present the eccentricities (top four panels) and the CDF of the eccentricity crossing orbit with respect to q (bottom panel) for each scenario.

As the mass of the planets increases, the gravitational influence of a stellar flyby causes greater destabilization within the planetary systems. This increase in planetary mass amplifies the perturbative effects, leading to more pronounced interactions among the planets. This increase in mass also results in a gradual growth in eccentricities for more distant flybys, reflected in the monotonic increase of the CDF (bottom panel of Figure 7). Each curve represents a distinct mass configuration, illustrating how the survival probability varies with the closest approach distance.

Furthermore, as planetary masses increase, their dynamical spacing decreases (measured in mutual Hill radii), impacting the system’s overall stability. For planets with masses approximately four times that of Jupiter, their Hill radius becomes comparable to that of Jupiter-mass planets in the 3:2-3:2 resonance. Consequently, the dynamical behavior of high-mass planets in the 2:1-2:1 resonance begins to resemble that of the 3:2-3:2 systems, shown in Figure 4. This is in agreement with the expectation that smaller dynamical spacings are an indicator of instability in compact configurations (Chambers et al. 1996).

As in the cases analyzed before, retrograde distant flybys increases the probability of the system retaining constant period ratios. In systems that do not cross the critical eccentricity threshold e_cross after the delayed integrations, the two- and three-planet resonant angles circulate, while the systems maintain constant period ratios without becoming unstable. Finally, in higher-mass systems, the difference between the immediate and delayed response to the passing star is significantly more pronounced, similar to what is observed in the 3:2-3:2 resonant chain.

Fig. 7 Four top panels show the delayed eccentricities for super-Jupiter systems starting in the 2:1-2:1 configuration. From top to bottom the set of planetary masses are (1.5, 2, 2.5) mJup, (2.25, 3, 3.75) mJup, (3, 4, 5) mJup and (3.75, 5, 6.25) mJup. Colors and markers follow the scheme used in Figure 4. The bottom panel shows the CDF for the distances of the flyby, each curve represents a different mass configuration.

3.3 The effect of the mass of the flyby star

In the previous section we showed that a 1 M_⊕ star approaching at distances of hundreds of astronomical units can break the resonant chains. Since most stellar encounters occur between unequal mass stars (Pfalzner 2013; Winter et al. 2018), we build on the fiducial case and vary the mass of the flying-by star, considering two different masses for the approaching star: M_fb = 0.5 and M_fb = 2 M_⊙. A similar approach is used for disks in Cuello et al. (2019), as noted in their Appendix D. For the planetary system, again we only focus on the case of the 2:1-2:1 resonant chain to analyze the overall impact.

The results are shown in Figure 8. As we established in the previous Sections, for the 2:1-2:1 resonance, the destabilization caused by the flyby is immediate, and for that reason we only show the results without the long-term additional integrations. Our simulations indicate that the larger the mass of the approaching star, the greater its destabilizing effect on the planets, in agreement with Malmberg et al. (2011). Moreover, for a more massive star, collisions with the central star and ejections from the system still occur even at larger relative distances, expressed as q/a₃. This suggests that the mass of the perturber (or the mass ratio between the host and intruding star) plays a crucial role in the dynamical evolution and stability of the resonant chains, should a flyby have occurred. However, based on stellar cluster dynamics and evolution, a 2 M_⊕ star approaching a solar mass star with planets is less likely to occur (Laughlin & Adams 1998; Pfalzner 2013).

It is important to emphasize, however, that while a close approach with a 2 solar mass star is not expected to occur frequently, the mass ratio between the stars in the simulations resembles that of GJ 876 encountering a 1 M_⊙ star. Conversely, when the flyby involves a star with M_fb = 0.5 M_⊕, the scenario is more similar to what HR 8799 might have experienced. This comparison allows us to place constraints on the flyby distance, assuming such an event occurred. In the left panels of Figure 8, we observe that if the minimum approach of the perturbing star exceeds 100 au from the planetary system, the system will most likely survive, which could explain the persistence of HR 8799-like systems, where such close encounters are less frequent. On the other hand, when a more massive star approaches the system, even distant flybys trigger collisions with the star and ejections, with fewer simulations resulting in constant resonant period ratios. A system such as GJ 876 could not have survived a flyby in a wide orbit configuration.

Fig. 8 Survival probability of period ratios of a three-planet system with masses (0.75, 1, 1.25) mJup in 2:1-2:1 resonant chains after a star with a mass of M* = 0.5 M⊙ (left) and M* = 2 M⊙ (right) flies by. In the bottom panels we show the cumulative distribution for the possible outcomes of the simulations.

4 Discussion

We investigate the orbital evolution of systems with long-period planets in resonant chains, focusing on the impact of a flyby star on their stability. Our analysis considers two initial configurations: 2:1-2:1 and 3:2-3:2 resonant chains.

The main findings can be summarized as follows:

1. Stellar flybys can disrupt wide-orbit resonant chains, either by ejecting planets immediately or through delayed interactions, potentially explaining the absence of such systems.

2. Jupiter-mass planets in 3:2-3:2 resonant chains are more easily disrupted by stellar flybys compared to 2:1-2:1 chains. As anticipated, the degree of disruption of the period-ratios increases when the star flies by the system in a retrograde orbit and with decreasing distance to the planetary system.

3. Flybys cause larger period ratio oscillations and force the resonant angles to transition out of the libration zone. Placing the outermost planet at 40 au, Jupiter-mass systems initially in a 2:1-2:1 resonant chain remain stable if the passing star remains farther than q = 125 au. On the contrary, systems in a 3:2-3:2 configuration already become unstable for distant flybys (q ~ 900 au) since they are at the brink of instability due to their more compact nature.

4. More massive planets experience disruptions across a wider range of flyby parameters, though increasing the planetary mass five times yields similar outcomes. Sub-Saturn mass planets maintain stable period ratios in both configurations, provided the flyby star’s mass remains constant.

5. For planetary masses up to that of Jupiter, the 2:1-2:1 chain responds instantaneously to flybys, whereas the 3:2-3:2 chain becomes unstable over time due to planet-planet interactions. The response of super-Jupiters to a stellar flyby in the 2:1-2:1 chain closely resembles the behavior of Jupiter-mass planets in the more compact 3:2-3:2 configuration.

6. The disruption of period-ratios commensurabilities is greater when a more massive star flies by the planetary system.

7. Jupiter-mass systems in wide 3:2-3:2 resonant chains are unlikely to survive stellar flybys, while lower-mass systems in wide orbits could remain stable.

Going back to the original question, it is possible that planetary systems originally formed in wide orbit resonant chains but were later disrupted by close encounters with passing stars. Such interactions could have destabilized the resonant configuration, resulting in the systems we observe today, where these chains no longer exist. If a system originally formed in a 3:2-3:2 resonance, disruptions seem more likely compared to those in a 2:1-2:1 configuration. In the case of HR 8799, the required q/a₃ value for disruption is sufficiently small to make this scenario unlikely.

One limitation of this study is the relatively narrow range of parameters explored. In our simulations, planetary systems orbit stars with fixed masses, unlike the variability observed in real systems. Varying the mass of the central star leads to differences in central potential, influencing how gravitationally bound the planets are to their host stars. These variations introduce complexities beyond the resonant dynamics. Future work could broaden the parameter space by considering a wider range of stellar and planetary masses, as well as orbital configurations. It is also important to note that while we propose this parametric study, the most expected outcomes are those resembling β = 45° or 135°. Stars passing with exact β values of 0°, 90°, or 180° are less likely to happen. Additionally, impulsive flybys interact with the planets for a shorter duration and may affect their orbits differently. However, such encounters are typically penetrating or highly catastrophic, likely leading to the complete disruption of the planetary system.

Finally, much remains to be explored concerning the influence of stellar flybys on the population of both detected and yet-to-be-detected exoplanets. In a resonant system, the associated angles are theoretically expected to librate. However, our simulations show that following a stellar flyby, these angles circulate instead, even as the period ratios remain constant. Thus, finding planets in such configurations – with stable period ratios but circulating angles – could potentially indicate that a past flyby event perturbed the planetary system.

Acknowledgements

The authors thank the anonymous reviewer for their report that helped improve our manuscript. C.C. acknowledges support from Agencia Nacional de Investigación y Desarrollo (ANID) through FONDECYT post-doctoral grant n°3230283. This project has received funding from the European Research Council (ERC) under the European Union Horizon Europe programme (grant agreement No. 101042275, project Stellar-MADE). C.P. acknowledges support from ANID BASAL project FB210003 and FONDECYT Regular grant 1210425. C.C. and C.P. gratefully acknowledge support from the John and A-Lan Reynolds Faculty Research Fund.

References

All Figures

	Fig. 1 Systems within 5% of resonant chains involving the 3:2 and 2:1 period ratio commensurabilities between adjacent planet pairs. Data have been taken from the Exoplanet Catalog (https://exoplanet.eu/catalog/). The size of the dots indicates planetary masses: large dots represent planets with masses of at least 1 mJup, while smaller dots correspond to less massive planets.
In the text

Fig. 2 Top: survival probability of period ratios (shading from lighter to darker) immediately after a star with a mass of Mfb = 1 M⊙ flies by a three-planet system with (0.75, 1, 1.25) Jupiter masses in 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains and leaves the surroundings of the system to a distance >1000 au, considering different inclinations and distances relative to the orbital plane of the planets. Bottom: cumulative distribution function (CDF) for the possible outcomes of the simulations: collisions with the star and ejections (red dashed line), survival of the period ratios (their resonant angles not necessarily librating, in solid purple), systems where period ratios do not remain commensurate (green), plotted against the pericentric distance of the flyby star relative to the outermost planet.

In the text

Fig. 3 Typical time evolution for the three planets in the 2:1-2:1 resonant chain after a close encounter with a flying-by star with mass 1 M⊙ and inclination β = 45°. The distance of minimum approach between the two stars is q = 125 au. In the top panels from left to right, we show the period ratio between the adjacent pairs, the eccentricities and inclinations, while in the bottom panels, the two-body resonant angles for each pair (left and middle plots) and the three-body resonant angle (in the bottom right panel).

In the text

Fig. 4 Delayed eccentricities (top row), inclinations (middle), and semi-major axis (bottom row) for all simulations in the 3:2-3:2 (left column) and 2:1-2:1 (right column) resonant chains. Colors indicate the pericentric distance q between the stars during the flyby event, with each color corresponding to a specific q, while the symbols represent the inclination β of the flyby star. The black dashed line indicates the stability limit given by the orbit-crossing value ecross.

In the text

	Fig. 5 Conditions leading to configurations within 5% of the 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains in the period-ratio plane (P2/P1, P3/P2). Relevant two- and three-planet MMRs are indicated by dashed gray curves. Same colors and markers as in Fig. 4.
In the text

	Fig. 6 Top panels: Period ratio survival for sub-Saturn mass planets, in 3:2-3:2 (left) and 2:1-2:1 (right) resonant chains. Bottom panels: delayed survival (after 105 orbits of the innermost planet) in the period-ratio plane. Colors and markers follow the scheme used in Figure 4.
In the text

Fig. 7 Four top panels show the delayed eccentricities for super-Jupiter systems starting in the 2:1-2:1 configuration. From top to bottom the set of planetary masses are (1.5, 2, 2.5) mJup, (2.25, 3, 3.75) mJup, (3, 4, 5) mJup and (3.75, 5, 6.25) mJup. Colors and markers follow the scheme used in Figure 4. The bottom panel shows the CDF for the distances of the flyby, each curve represents a different mass configuration.

In the text

	Fig. 8 Survival probability of period ratios of a three-planet system with masses (0.75, 1, 1.25) mJup in 2:1-2:1 resonant chains after a star with a mass of M* = 0.5 M⊙ (left) and M* = 2 M⊙ (right) flies by. In the bottom panels we show the cumulative distribution for the possible outcomes of the simulations.
In the text