On TTVs with rebound

[douglasalvesastro12]

Hello @hannorein thank you for making rebound available to us, it's an amazing code!

As an exercise for a larger project I was trying to make some TTVs following the python examples and I can't reach conclusive results. Hopefully you can help me with the following:

ex1: For a 3 body system, M0 = 1 Msolar, m1 = 1e-3 Ms(~Jupiter) at P1 = 12 days and m2 = 3*1e-6 Ms (~Earth) at P2 =16 days, we would have a roughly sinusoidal TTV shape with amplitude ~4min. Both planets are in a coplanar orbit with zero eccentricities.

What I get is something completely different (image attached below).

• The system is in 1st order MMR j:j+1, j=3 and I believe I have the units consistent.
  Questions:
  (0) Why I am not getting the sinusoidal shape with amp. ~ 4min?

(1) What are reasonable numbers for these lines: System integration sim.integrate(sim.t+1e-4), finding the central time while t_new-t_old>1e-7 and getting out of central time sim.integrate(sim.t+1e-3). I used 1e-4, 1e-7 and 1e-3 which is in years and represents respectively ~<1% of period, left as is and ~3% of period. These numbers if not chosen correctly seems to heavily influence both the shape and amplitudes of TTVs.

(1) I noticed that changing omega (argument of pericenter) the TTV shapes change. Should it not be the same, specially for circular orbits? Circular orbits don't have apogeo/perigeo hence omega is not actually defined for circular orbits. (right?). However, even for the case of 1 planet in inner circular orbits and an outside perturber in a non-circular orbits, omega should not make a huge difference in the TTV shape of the inner because in the long term there will be a precession of omega in a somewhat periodic manner.

The snippet I used from your examples

fig = plt.figure(figsize=(10,5))
ax = plt.subplot(111)

        
sim = rebound.Simulation()
sim.G = 39.478 #AU^3 yr^-2 Ms^-1
sim.add(m=1.)
sim.add(m=1e-3, P=12./(365.25),e=0.0)
sim.add(m=3*1e-6, P=16./(365.25), e=0.0)
sim.move_to_com()

N=50
transittimes = np.zeros(N)
p = sim.particles
i = 0
while i<N:
    y_old = p[1].y - p[0].y 
    t_old = sim.t
    sim.integrate(sim.t+1e-4) # check for transits every 0.5 time units. Note that 0.5 is shorter than one orbit
    t_new = sim.t
    if y_old*(p[1].y-p[0].y)<0. and p[1].x-p[0].x>0.:   # sign changed (y_old*y<0), planet in front of star (x>0)
        while t_new-t_old>1e-7:   # bisect until prec of 1e-5 reached
            if y_old*(p[1].y-p[0].y)<0.:
                t_new = sim.t
            else:
                t_old = sim.t
            sim.integrate( (t_new+t_old)/2.)
        transittimes[i] = sim.t
        i += 1
        sim.integrate(sim.t+1e-3)#0.1)       # integrate 0.05 to be past the transit

A = np.vstack([np.ones(N), range(N)]).T
c, m = np.linalg.lstsq(A, transittimes, rcond=-1)[0]

plt.plot(range(N), (transittimes-m*np.array(range(N))-c)*(60.*24.*365.25), '.-')
ax.set_xlim([0,N])
ax.set_xlabel("Transit number")
ax.set_ylabel("TTV [min]")

[hannorein]

PSA: The bisection method used in the REBOUND example is a bit crude. I'm currently working with a summer student to hopefully come up with a slightly more elegant way to do this. (But there's nothing wrong with it.)

The numbers that you have chose seem to be reasonable for your case.

• The first one (1e-4) needs to be smaller than the orbital period otherwise it might miss a transit.
• The second one (1e-7) determines the accuracy. The smaller the more accurate. It should be smaller than the accuracy that you want.
• The third one (1e-3) makes sure you don't detect a transit twice. Anything longer than the transit duration and shorter than an orbital period is fine.

Having zero eccentricity is a very special case when it comes to orbital elements. They way orbital elements are implemented in REBOUND ensures that they are "smooth" when the eccentricity goes to zero. So the omega will matter even in that case.

Are you sure the planets are in a MMR? The period ratio is an integer ratio, but the eccentricities are zero and one planet is much more massive than the other. I suspect the plot that you made is the correct TTV signal. Maybe you're already doing this, but I'd start with trying to reproduce a well known TTV system. (If that's what you're doing, can you share which system this is?)

[douglasalvesastro12]

Thank you for your reply @hannorein, good to know the numbers I used were ok.
Regarding your statement: "They way orbital elements are implemented in REBOUND ensures that they are "smooth" when the eccentricity goes to zero. So the omega will matter even in that case."
How may I choose omega? The code I posted gives that TTV shape and amplitude within the 1min level, however, if I set omega=90 here sim.add(m=1e-3, P=12./(365.25),e=0.0) I get a different shape with amplitudes ~60min hence choosing omega wrong severely impacts the TTV shape. I notice that changing the precision from 1e-7 to 1e-10 with omega=90 for planet 1 also changes the TTV shape. Do you know why these things happen?

"Are you sure the planets are in a MMR? The period ratio is an integer ratio, but the eccentricities are zero and one planet is much more massive than the other."
It is supposed to be in a 1st order MMR 3:4. P2/P1 = 4/3, but as this is an example (not a real case) I should try to reproduce a real example, you're absolutely right. I will check on a real case (Kepler-9, for instance) and let you know. Really thanks @hannorein

One more question/comment: REBOUND does not seem to care about stellar-to-planet ratio when transits occur, you check roughly when the planet crosses the y-axis. This puts an uncertainty in the TTVs computed with REBOUND which would be in the same level of the transit duration, right? Is there a way to go around it or would you guys implement that? Perhaps asking for these (approximated) radius and then checking for the central time within the regime of |y| < Rstar ?

[LefterisDimitri]

Hello @hannorein and @douglasalvesastro12 ,

I have a question about the number of transit points N. Does this number has a connection with the period of the planet that we study?

Thank you

[douglasalvesastro12]

Hi Dimitriou, n is usually used to represent the transit number which relates to the planet period and central transits as Tn = T0 + n*P, this way you can get the central time of your transits (assuming zero uncertainty in the first mid transit and planet period P) Cheers

…

On Tue, Jun 27, 2023, 12:28 Dimitriou Eleftherios ***@***.***> wrote: Hello @hannorein <https://github.com/hannorein> and @douglasalvesastro12 <https://github.com/douglasalvesastro12> , I have a question about the number of transit points N. Does this number has a connection with the period of the planet that we study? Thank you — Reply to this email directly, view it on GitHub <#551 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AI5VVSSZOBFUHKB65K6SDEDXNK7XVANCNFSM4654MYTQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[LefterisDimitri]

Thank you @douglasalvesastro12