Create comets from particles in mean motion resonance with an eccentric planet

[tastuber]

Environment

• REBOUND Version: 4.4.1
• API interface: Python
• Operating System (including version): MacOS 14.6.1 or CentOS Linux 7 (on HPC cluster)

Physics
I want to simulate a planetary system with planetesimals in mean-motion resonance with an eccentric planet. The secular influence of the planet causes the eccentricity of the planetesimals oscillate. When the eccentricity is high, the orbits of the planetesimals cross the orbit of the planet and can get scattered onto cometary orbits during a close encounter.
Over the timescale of ~1Gyr I want to monitor cometary activity in such a system.

To begin with I want to reproduce the results from Faramaz et al. 2017 who investigated this problem numerically with the code SWIFT-RMS. That code is fully symplectic with an adaptive time step.

I chose a particular simulation to reproduce:

• The central star is solar like.
• The planet of mass $0.1 \textrm{M}_\textrm{Jup}$ is on an orbit with semi-major axis 100au and eccentricity 0.1
• There are 5000 planetesimals.
• The planetesimals have
  • Semi-major axes from 54.04au to 54.54au
  • Eccentricities from 0.0 to 0.05
  • Inclinations from -3deg to 3deg
  • The other orbital elements are picked from a uniform random distribution.

This distribution of planetesimals encloses the 5:2 mean motion resonance with the planet. It is not meant to reproduce a particular system, but to explore the above described mechanism to produce comets.

Goal
Investigate the production and lifetime of comets produced by the above described mechanism over the timescale of Gyrs.
Final measure:

• number of cometary events per year
• cometary lifetime.

Progress
Due to the nature of the problem with long integration times (~ Gyr) paired with the need to accurately compute close encounters I chose the hybrid integrator Mercurius.

Using the webinterface I checked that the simulations are set up as I expected.

Over a time of ~20Myr the planetesimals in the mean-motion resonance are excited up to eccentricities of ~0.8. From rough analytical estimations this should be correct.

I analyzed the simulation archive to detect active comets. To do that I assumed a particle to be a comet if its periastron is less than 3au.
The results I get are not like in Faramaz et al. 2017 where comets start to appear after ~500Myr (their Fig. 5).
When I use a integration timestep of 5% of the smallest orbital period, I get comets right after the planetesimals got excited, i.e., around ~20 - 25Myr. If I choose a timestep of 1% of the smallest orbital period I do not get any comets within 1Gyr.

Question
I think the scattering behavior of the planetesimals is strange. I can elaborate more on that if necessary.
But my first question is: did I set up the simulation correctly?

import rebound
import astropy.units as u
import numpy as np

# Set up random number generator and seed it.
rng = np.random.default_rng(2024)

# Set final integration time and snapshot interval.
t_final = 1e9 #  [yr]
dt_sim_archive = 1e6 #  [yr]

# Use a fraction of the shortest orbit as the integration time step.
frac_shortest_orbit = 0.05

# Star properties (central mass).
m_star = 1.0 #  [Msun]

# Planet properties.
m_planet = 0.1 * u.Mjup.to(u.Msun) #  [Msun]
a_planet = 100.0 #  [AU]
e_planet = 0.1

# Planetesimal belt properties; "tp" stands for test particles.
N_tp = 5000
a_tp_min = 54.04 #  semi-major axis
a_tp_max = 54.54
e_tp_min = 0.0 #  eccentricity
e_tp_max = 0.05
inc_tp_min = -3.0 * (np.pi/180) #  inclination [rad]
inc_tp_max = 3.0 * (np.pi/180)

##### Set up rebound simulation. #####

sim = rebound.Simulation()
sim.units = ('yr', 'AU', 'Msun')
sim.integrator = 'mercurius' #  hybrid integrator
sim.ri_ias15.min_dt = 0.0
sim.testparticle_type = 0 #  massless non-interacting testparticles

# Set up simulation archive.
sim.save_to_file(
    "testarchive",
    interval=dt_sim_archive, delete_file=True
)

# Add active bodies, i.e., the star and the planet
star = rebound.Particle(m=m_star)
sim.add(star)
planet = rebound.Particle(
    simulation=sim,
    m=m_planet,
    a=a_planet,
    e=e_planet,
    f=0.0
    )
sim.add(planet)
sim.N_active = sim.N #  sim.N = 2

# Add planetesimal belt.
# Initialize the belt by drawing uniformly distributed random numbers from the
# set intervals.
while sim.N < (sim.N_active + N_tp):
    a = rng.uniform(a_tp_min, a_tp_max)
    e = rng.uniform(e_tp_min, e_tp_max)
    inc = rng.uniform(inc_tp_min,
                      inc_tp_max)
    Omega = rng.uniform(-np.pi, np.pi)
    omega = rng.uniform(-np.pi, np.pi)
    M = rng.uniform(-np.pi, np.pi) #  mean anomaly
    p = rebound.Particle(
        simulation=sim,
        a=a, e=e, inc=inc, Omega=Omega, omega=omega, M=M
    )
    sim.add(p)

# Compute the integration time step.
orbits = sim.orbits()
shortest_orbit = np.inf
for orbit in orbits:
    if orbit.P < shortest_orbit:
        shortest_orbit = orbit.P
timestep = frac_shortest_orbit * shortest_orbit
sim.dt = timestep

# Move to center of mass frame.
sim.move_to_com()

# Compute total energy before the simulation.
E0 = sim.energy()

# Perform the integration.
sim.integrate(t_final, exact_finish_time=0)

print("Integration finished.")

# Print out relative energy error after integration.
E1 = sim.energy()
print(f"Relative energy error: {(E0-E1)/E0:.2e}")

[hannorein]

Hi tastuber (what's your real name?),

Your setup looks alright. But in general this can be a difficult simulation to get right.

• Try out @tigerchenlu98's new TRACE algorithm. This should be able to handle this setup better.
• Check if any of the test particles come close to either the star or the planet. If so, you might need to include the physical size of these objects and think about how to handle collisions.

[hannorein]

(I guess you're Thomas.. sorry... didn't make the connection immediately)

[tastuber]

Yes, I am Thomas!

Thank you, I will try out the TRACE algorithm.
It also used the high-order IAS15 integrator for close encounters. Do you expect any impact on the long time accuracy due to the non-symplectic nature of this integrator (or mercurius)?
Are there reasons arguing against a fully symplectic integrator for my investigation?

I would have guessed that collisions are extremely rare and would not change the results, which are of statistical nature anyway.
Do you expect any significant impact of collisions?

There is another thing I do not consider at the moment:
After so and so many periastron passages a comet is gone (either completely evaporated or broken up). One would have to compare the typical survival timescale of evaporation against the the dynamical survival timescale of the cometary orbit itself.
I expect that the impact of this would be larger than considering collisions.
For now I am primarily interested in whether this setup produces comets and if yes, when.

[hannorein]

IAS15 is as accurate as one any integrator can get on a computer with double floating point precision. Therefore it is also as symplectic as any symplectic integrator. So no, this is not an issue.

If physical collisions are rare, then, no, they should not be important.

One of my students, @dangcpham has recently run similar simulations. As in his case, your planet and star are not affected by the 5000 particles. This can be uses to speed up these simulations. But this is not implemented in the default version of REBOUND. I'm chaining him in here in case he wants to elaborate.

[tastuber]

Do you have any suggestions about the time step?
I saw a striking difference in scattering frequency when I changed from 5% of the shortest orbital period to 1%.
During some simulations (different planet masses) I got the notification that the time step is larger than at least one orbital period.
I can only explain this to me that some particle got scattered into a very short orbit while still being outside the x Hill radii.

The results I am comparing mine with were computed using a fully symplectic integrator with an adaptive timestep. If a particle got closer to the planet/star than x Hill radii, the time step was reduced.
Would such a scheme produce results of comparable precision compared to TRACE or Mercurius?

[hannorein]

Have you tested if you get different results using TRACE?

The 5% timestep should be fine. And yes, you should be able to reach the same precision.

You can always run the simulation with IAS15 to get the most accurate results. It'll be slower, but you can for example reduce the number of particles.

[dangcpham]

Hi Thomas,

This might be a lot so I'm also happy to meet on Zoom if you want to discuss this further!

1. IAS15 is as accurate as computers can allow it, so although it is non-symplectic it shouldn't matter.
2. I don't see where you're handling collisions with the star or the planet. One solution is using a C heartbeat function at every timestep and check if the comet's distance is within the star's or planet Hill radius.
3. I have run a similar setup before, but since particles massless and non-interacting, I just run multiple star-planet-comet simulations and put all the results together at the end. I find that simulations in this case run faster as they are adapting timesteps more appropriately.
4. Your notification that the time step is larger than at least one orbital period is probably due to IAS15 proposing a large timestep at some point, but then the comet decreased in semi-major axis and is now interior to the previously most interior comet. Then Mercurius switched over to WHFast but didn't update with the innermost orbit. Using TRACE should fix that. Alternatively, you can just run purely on IAS15 and that should never be a problem.
5. If your interests is only star-planet(s)-comet, then the numerical method in my recent paper might help you run this much faster: https://arxiv.org/pdf/2404.07160. See Figure 8, and appendices A, B, C where I showed the precisely your case: the Kuiper Belt interacting with a planet. It runs a lot faster and still reproduces resonances and planet-comet interactions correctly.
6. I wouldn't worry about the chemical evolution. Planets can be quite effective at strongly scattering comets. Previously works like (https://ui.adsabs.harvard.edu/abs/1997Icar..127...13L/abstract, https://ui.adsabs.harvard.edu/abs/2023AAS...24113602K/abstract) typically care about comet fading as they relate to comet-planet interactions, not from volatiles depletion. In the Solar System, I think only the Near-Earth Asteroids (3-5 au) are comets depleted of volatiles. These have much closer orbits than your 30-50 au comets and thus, are constantly sublimating materials in contrast to yours.
7. If you want to show (6) quantitatively, you can calculate the chemical lifetime for an eccentric orbit from Table 5.1 and Chapter 5 in (https://www.issibern.ch/PDF-Files/SR-004.pdf). Then, compare that with the dynamical lifetime from your N-body simulation.

[tastuber]

Hi Dang,

1. I will run simulations for comparison with IAS15.
2. So far I do not consider collisions at all.
   This is mainly because I want to compare my results with those of a specific setup (Faramaz et al. 2017) where collisions are neglected as well.
3. Did you combine the Simulation Archives to put the results together? Is there such a functionality?
4. My setup is slighlty different: the planet orbit is outside the planetesimal orbit. But I will try TRACE and IAS15.
5. I briefly skimmed through your paper and the setup sounds very neat. For now it would probably be overkill in my case. FIrst I have to get stable results.
6. & 7. I should have defined my usage of the term 'comet'. I was considering a planetesimal as a comet only when its periastron was close enough to start rapid sublimation. I usually assume 3au aroung a solar like star; here they note 2au Marboef et al. 2016.

I will run simulations with TRACE, IAS15, and MERCURIUS for comparison and report back with the results.
Thank you for the Zoom offer, I will keep it in mind!

[dangcpham]

3. I kept track of the comets' evolution and then process all of them together afterward.

4. The method actually is pretty easy to implement, so let me know if you want to use it later!

6&7. In the case of rapid sublimation near periastron, the sublimation rate is actually relatively straightforward and you can roughly calculate how much is lost over one periastron passage (see Equations 3.6 and 3.7 in https://link.springer.com/book/10.1007/1-4020-3495-4). The 2-3 au number is fascinating. I'm currently finishing a comet chemical evolution paper myself and tracked it down. It turns out to be ~1 au for water ice, which is the majority of cometary volatiles. You can read a bit more about that in the 2005 book I linked earlier, or the original Delsemme chapter in this 1981 Comets book (https://archive.org/details/comets0000unse). Good luck!

[tastuber]

Do you recommend to change any TRACE settings from the default?

[hannorein]

I think the defaults should work fine.

[dmhernan]

@tastuber , not that it's too relevant to the discussion, but SWIFT-RMVS works by switching the Hamiltonian splitting every so often, as does TRACE. And as we showed in appendix of Lu et al. (2024), this means neither TRACE nor RMVS are actually fully symplectic. And in practice, this can degrade your error over time (see, e.g., Fig. 7 here https://arxiv.org/pdf/2301.06253). But TRACE is nearly time-reversible, and as a result doesn't degrade in error.

Probably not information you needed, but all to say I'd only expect improvements from using TRACE in your case.

[tastuber]

Thank you for the information!
It might help interpreting new and past results.

My TRACE and IAS15 simulations are still running.
I guess is is expected that TRACE takes longer than Mercurius.