migration to a target value

[nchoksi4456]

Hi Hanno et al. I hope you have been well. Is there an easy way to implement migration to some target value? I'm trying to implement something of the form da/dt = -(a - a_target)/tau_a.

Thanks!
Nick

[hannorein]

Maybe something like this?

https://github.com/dtamayo/reboundx/blob/main/ipython_examples/ExponentialMigration.ipynb

[nchoksi4456]

Thanks. Yea, something like this. But I've attached an example that showcases an issue I'm having. The example should run out-of-the-box.

Basically, I'm trying to set up a three planet system with m1 = 3e-5, m2 = 0, m3 = 3e-5. The goal is to establish a chain of near-resonances. In this goal configuration, m2 is 1% away from 3:2 commensurability with planet 1, and 1% away from 3:2 commensurability with planet 3. I'm trying to use ExponentialMigration to drive m2 to the orbital period that satisfies this configuration. I think it is getting close. But, not close enough. Since we're working near resonances, tiny (percent-level) deviations in orbital period from my desired value make a big difference. See attached code and plot, which shows that m2 is driven a half-percent or so shy of the target orbital period.

I've confirmed that when I set the other bodies to be test particles, ExponentialMigration drives the test particle to exactly the right target orbital period. Not sure what (if anything) can be done about this. Thoughts welcome!

import rebound 
import reboundx
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import time as ti
def setup_and_run_sim( fname, args, density = 1e4, useP = True):
    sim = rebound.Simulation()
    sim.integrator = "mercurius" # Integration method. Wisdom-Holman is good for simple, periodic orbits. Not good for close encounters.



    sim.G = 1
    sim.add(m = 1)

    Np = args["Np"]
    # NC: Assign planet properties using the args dictionary
    for idx in range(1, Np+1):
        sm = "m%d" % idx
        if useP:
            sP = "P%d" % idx
        else:
            sa = "a%d" % idx
        se = "e%d" % idx
        so = "omega%d" % idx
        sf = "f%d" % idx
        if useP:
            sim.add(m = args[sm], P = args[sP], e = args[se],  omega = args[so], f = args[sf])
            #print("Added particle of mass = %e , period P  = %.3f, eccentricity e %.3f " % (sim.particles[-1].m, sim.particles[-1].P, sim.particles[-1].e))
        else:
            sim.add(m = args[sm], a = args[sa], e = args[se],  omega = args[so], f = args[sf])
            #print("Added particle of mass = %e , semimajor axis a  = %.3f, eccentricity e %.3f " % (sim.particles[-1].m, sim.particles[-1].a, sim.particles[-1].e))

    sim.dt = 0.05*sim.particles[1].P
    
    rebx = reboundx.Extras(sim)
    mof = rebx.load_force("modify_orbits_forces")
    rebx.add_force(mof)
    
    mod_effect = rebx.load_force("exponential_migration")  # Add the migration force
    rebx.add_force(mod_effect)  # Add the migration force
    
    

    # NC: Assign migration and eccentricity damping times
    ps = sim.particles
    for idx in range(1,Np+1):
        ste = "te%d" % idx
        sta = "ta%d" % idx

        if(abs(args[sta]) < 1e15/sim.particles[1].n):
            print('adding for idx = %d' % idx)
            ps[idx].params["em_aini"] = ps[idx].a
            ps[idx].params["em_afin"] = args['a_target']
            ps[idx].params["em_tau_a"]  = args[sta]

        if(abs(args[ste]) < 1e15/sim.particles[1].n):
            ps[idx].params["tau_e"] = -args[ste]


    sim.move_to_com()
    start = ti.time()
    
    sim.save_to_file(fname, interval=  100, delete_file=True)
    
    
    # NC: Integrate until the middle test particle reaches the desired semimajor axis
    achieved_target = False
    while achieved_target == False:
        sim.integrate(sim.t + 1000*sim.dt, exact_finish_time = 0 )
        if sim.t > (10*abs(sim.particles[2].params["em_tau_a"]) ):
            break
    
    
    rebx.remove_force(mof)
    rebx.remove_force(mod_effect)

    
    
j = 3

# NC: Three planet system with m1 = 3e-5, m2 = 0, m3 = 3e-5
# NC: Goal is to establish a 3 planet, near-resonant chain
# NC: In this goal configuration, m2 is 1% away from 3:2 commensurability with planet 1, and 1% away from 3:2 commensurability with planet 3


P1 = 2*np.pi
D12_target = 0.01
P2_target = P1*(j/(j-1))*(1 + D12_target)
a2_target = (2*np.pi/P2_target)**(-2./3)  # NC: Desired semimajor axis
D23_target = 0.01
P3 = (j/(j-1))*P2_target*(1 + D23_target)


P2_initial = 1.1*P2_target  # NC: Initial orbital period != desired orbital period P2_target. Goal is to have migration drive the test particle to P2_target

te = 1e5

args = { "Np":3, "te1":te, "te2":te, "te3":te, "ta1": 1e20, "ta2":1e4, "ta3": 1e20 , "m1":3e-5, "m2":0, "m3":3e-5, "P1":P1, "P2": P2_initial , "P3":P3, "e1":0, "e2":0, "e3":0, "omega1":0.8, "omega2": 0.8, "omega3":4, "f1":0, "f2":3.9, "f3":1}
args['a_target'] = a2_target


setup_and_run_sim("test", args)
sa = rebound.Simulationarchive("test")

tlist = np.array([sim.t for sim in sa])
P2list = np.array([sim.particles[2].P for sim in sa])
plt.plot(tlist, P2list/P2_target)
plt.ylim(0.98, 1.02)

plt.ylabel("P2/target P2", fontsize=19)
plt.xlabel("time (1/n)", fontsize=19)

[hannorein]

This seems more like a physics question. In the original Lee and Peale papers, the migration doesn't stop when the planets reach commensurability. It's the resonance itself which locks the period ratio to be a constant. And once you look at the 0.1% level like in your plot, then all kind of other things become important. For example, I'm not sure if heliocentric or Jacobi coordinates are the right choice here. Looking at resonant angles might in general be a better diagnostic for the resonance than a period ratio.

[dtamayo]

Hi, like Hanno said this is really tricky. We had to worry about this for our 2017 Trappist paper too. Our approach was to just run migration simulations and then use the scale invariance of gravity to rescale everything to the inner period observed and keep all the period ratios / resonance relations. You can find code here if that sounds like it would be helpful

https://github.com/dtamayo/trappist