Confusing function name in "type_I_migration.c"

[miguelcuenro]

Hi,

in the file "type_I_migration.c", that implements type I migration, what is called the semi-major axis damping timescale is probably meant to be the migration timescale, which isn't the same. This can be a source of confusion and it might be a good idea to use a different naming convention (e.g. migration timescale, as in Cresswell&Nelson2008).

[hannorein]

What's the difference?

[miguelcuenro]

Hello again!

First of all thanks for the fast reply!!

The migration timescale is equal to the angular momentum divided by the torque (\tau_mig = L / \dot{L}), while the semi-major damping timescale is equal to the semi-major axis divided by its time derivative (\tau_a = a / \dot{a}), which results in:

1/tau_a = (1/(tau_mig/2) + C(e)/(tau_e/2)), C(e) = e^2/(1-e^2)

In other words, the migration timescale is the damping timescale of the orbital angular momentum, while the semi-major axis damping timescale is the damping timescale of the semi-major axis, and the two are not equal to each other (not even in the circular case, there is a factor 2).

Therefore, it could be convenient to use a non-conflicting name, such as "rebx_calculate_migration_timescale", rather than "rebx_calculate_semi_major_axis_damping_timescale".

[hannorein]

Thanks! I can see that the difference between angular moment and semi-major axis damping timescale can be important. But I'm less sure about why one of them deserves the name migration timescale but not the other. I guess you consider Cresswell & Nelson (2008) to be the authority for this?

[miguelcuenro]

Yes, I did take Cresswell & Nelson (2008) as my main source, since the migration timescale, defined in equation 13, is the one implemented by "type_I_migration.c".

[dtamayo]

Thanks so much for helping clarify the documentation to make it easier for everyone else!

That sounds reasonable to me. I've been traveling today and can't go through it carefully right now.

Are you proposing to rename

rebx_calculate_semi_major_axis_damping_timescale

Because it's currently calculating the L/Ldot 'migration' timescale

Or to rename the currently implemented

rebx_calculate_damping_timescale

Does that sound reasonable to you @acpetit ?

[GabrielePichierri]

Hi all, if I may jump in the conversation, Miguel and I want to implement other migration prescriptions into reboundx so we just want to make sure we're on the same page with what's there already. I would say rebx_calculate_damping_timescale can stay as it is (it enters in all damping timescale formulas). The issue is more with rebx_calculate_semi_major_axis_damping_timescale, which doesn't calculate tau_a but rather tau_mig=L/ \dot{L} (Miguel tested this running a simulation). And since migration is normally associated with a torque (dot{L}) provided by the disk, it would then make sense to call it something like rebx_calculate_migration_timescale, plus as you say it keeps in line with Cresswell & Nelson.

[hannorein]

Thanks! I can see what the issue is now!

[GabrielePichierri]

Great! Antoine is on vacation at the moment but he says he agrees.

I think other than in modify_orbits_direct, which clearly modifies the orbital elements by design, all other (more realistic) migration prescriptions should refer to the angular momentum damping timescale (tau_mig) rather than tau_a. This is because all the literature on the subject of planet migration focuses on the torque, so people (like us) will want to implement torque prescriptions.

[dtamayo]

Thanks all, I understand now. So the problem is that we need to write a function that calculates the x,y,z components of the acceleration on each particle. Right now type_I is reusing the code from modify_orbits_forces, which implements the equations of Papaloizou & Larwood to apply a force opposite the overall velocity, and another opposite the radial component of the velocity. The coefficients are chosen so that the semi major axis and eccentricity damp exponentially at the rate specified by tau_a and tau_e.

Clearly the same can be done by matching an energy and angular momentum loss. Is that something you (Miguel or Gabriele) could easily code up in src/type_I_migration.c from Cresswell & Nelson or some other source? I'm happy to do the renamings at the appropriate points in the code.

[GabrielePichierri]

I agree that modify_orbits_forces.c does what the name suggests.(*) I'm not sure however the code type_I_migration.c reuses the code from modify_orbits_forces.c: It contains itself the implementation of Papaloizou & Larwood eqs (37-38); in particular

    if (invtau_a != 0.0){
        a.x = -dvx*(invtau_a);
        a.y = -dvy*(invtau_a);
        a.z = -dvz*(invtau_a);
    }

implements eq (37), only what is called invtau_a is actually what Papaloizou & Larwood call 1/t_m, which they define as \dot(angular momentum)/(angular momentum).

(*) Actually this might be strictly true in the circular case or when there's no e-damping, because right now in modify_orbits_forces.c it's written:

    a.x =  dvx*invtau_a/(2.);
    a.y =  dvy*invtau_a/(2.);
    a.z =  dvz*invtau_a/(2.);

In other words it's defining the 1/t_m of Papaloizou & Larwood as 1/t_m=invtau_a/(2.), but this is true only when e=0 or 1/tau_e=0.

[GabrielePichierri]

(P.S.

(*) Actually this might be strictly true in the circular case or when there's no e-damping

I checked this and it seems to be the case. However you state something of the sort in the documentation so the user should already be aware of it.)

[dtamayo]

I see. Is the problem then just in the names used in the source code? The user never sets a 'tau_a' or a 'tau_mig' for the type I implementation right? But I agree it avoids future problems for future effects with different torque prescriptions.

It sounds like Antoine and you are in agreement, so perhaps it would be most efficient for you to make the changes you think are needed and send a pull request? If you let me know what you think needs to change, I'm also happy to make the changes, but it would have to be next Tuesday as I'm traveling early tomorrow.

Thank you for helping improve this!

[GabrielePichierri]

Yes I think the user who just wants to use to code as it is doesn't even see the tau_a' or 'tau_mig', so in the end I can just make the name change and nobody should be affected. I'll do the pull request.

Cheers!

[dtamayo]

This was updated in Gabriele's pull request now in 3.7.2. Closing