Merging Issues with a Large Number of Particles

[PhamNguyen18]

Hello,

I am interested in using REBOUND to model the Moon forming accretion disk but have encountered an issue with particle mergers. As a test, I would like to replicate one of the simulations from Ida 1997. I am using the planetesimal disk migration example file as a template for my simulation and parameters from run #4 in Ida’s paper. Below is some of the code I am using to initialize the disk (sorry for the somewhat long length):

    // With in the main function
    r->integrator    = REB_INTEGRATOR_MERCURIUS;
    r-> G = 6.67430e-11; //KMS

    // Collisions
    r->collision = REB_COLLISION_DIRECT;
    r->collision_resolve = reb_collision_resolve_merge_pass_through;
    r->track_energy_offset = 1;
    r->collision_resolve_keep_sorted = 1;

    // Constants
    double m_earth = 5.9722e24; // kg
    double r_earth = 6.371e6; // m
    double m_lunar = m_earth * 0.0123;
    double aroche = 1.85518e7; // Roche radius
    double t_kepler = 25200; // Keplerian orbital period at roche (7 h)

    srand(12);
    r->dt = 0.002*t_kepler;

    // Earth 
    struct reb_particle earth = {0};
    earth.m         = m_earth;
    earth.r        = r_earth;    
    reb_add(r, earth);

    // Planetesimal disk parameters
    double total_disk_mass = 2.44*m_lunar;
    int N_planetesimals = 50;
    double planetesimal_mass = total_disk_mass/N_planetesimals; 
    double powerlaw = 0.25;
    double amin = 0.35*aroche, amax = aroche*0.95; 

    // Generate Planetesimal Disk
    while(r->N<N_planetesimals + 1){
        struct reb_particle pt = {0};
        double dist = 0;
        // Make sure no particles are initialized inside the Earth!
        while (dist < r_earth*r_earth) {
          double  a = reb_random_powerlaw(amin,amax,powerlaw);
          double e    = reb_random_rayleigh(0.2); 
          double inc  = reb_random_rayleigh(0.2);
          double Omega = reb_random_uniform(0,2.*M_PI);
          double apsis = reb_random_uniform(0,2.*M_PI);
          double phi     = reb_random_uniform(0,2.*M_PI);
          pt = reb_tools_orbit_to_particle(r->G, earth, planetesimal_mass, a, e, inc, Omega, apsis, phi);
          dist = pt.x*pt.x + pt.y*pt.y + pt.z*pt.z;
        }
        pt.r         = particle_r(planetesimal_mass, m_earth, aroche);
        reb_add(r, pt);
    }

    r->N_active = r->N;

  // Outside main function
  // Used to calculate plantesimal radius
  double particle_r(double m, double m_earth, double roche_radius) {
    return (1.0/2.456) * pow(m/m_earth, 1.0/3.0) * roche_radius;
  }

To keep things simple, I have kept the number of particles small for now (run #4 from Ida uses 2000 particles) and of uniform mass. Besides the planetesimal disk parameters, the main modification I have introduced is allowing for gravitational forces between all particles in the disk and the Earth (setting N_active = N). My issue is that after a brief merging period the simulation returns “nan” for the relative energy difference and for the particle positions as well! I am new to using REBOUND and was hoping to get some advice on how to address this problem. Thank you!

[hannorein]

I tried running the example. It looks like all the particles just merge into one after a short time. I don't get any "nan"s. When I reduce the particle radii, nothing collides anymore.

• Are you using the latest version of REBOUND?
• Are you sure you setup the radii correctly? The radii are very large (see screenshot).

[PhamNguyen18]

Hello Dr. Rein,

I updated REBOUND just in case. I think you are right that the radii are too large. The function I used to calculate the radii is from Ida but I think using only 50 particles leads to planetesimals that are too large in radii. In the actual paper they range in mass from 10^-5 to 10^-2 lunar masses and follow a power-law size distribution. I’ll decrease the particle size and see what results I obtain with a larger number of particles as well.

[hannorein]

I'll close this for now. But if you the problem persists, feel free to reopen the issue.