Reproducing real TTVs with ttvfast-python

[douglasalvesastro12]

Hello,

I am trying to reproduce the TTVs from (Holman, M., J, et al 2010): Kepler-9 - A system of Multiple Planets Transiting a Sun-like Star, Confirmed bu Timing Variations. A google image of the TTVs can be seen hopefully here (top figure you see Kepler9-b in blue and kepler9-c in red) https://www.google.com/imgres?imgurl=https%3A%2F%2Fscience.sciencemag.org%2Fcontent%2Fsci%2F330%2F6000%2F51%2FF3.large.jpg&imgrefurl=https%3A%2F%2Fscience.sciencemag.org%2Fcontent%2F330%2F6000%2F51&tbnid=mZIzIzB4zCNX6M&vet=12ahUKEwjwtZrwzbrxAhUXkZUCHbPLAkoQMygDegUIARCWAQ..i&docid=Xlb5WQ8q-INRpM&w=755&h=1280&q=Kepler-9%3A%20A%20System%20of%20Multiple%20Planets%20Transiting%20a%20Sun-Like%20Star%2C%20Confirmed%20by%20Timing%20Variations&ved=2ahUKEwjwtZrwzbrxAhUXkZUCHbPLAkoQMygDegUIARCWAQ. The idea is to use TTVfast-python for a different target for a paper, however I am struggling a bit to understand the code behavior. The piece of code I am using follows:

import ttvfast
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

# gravity = 0.000295994511                        # AU^3/day^2/M_sun
stellar_mass = 1.0#0.95573417954                    # M_sun

#planet parameters come mostly from exofop and posterior dist. from previous modeling. except for longnode, argument, mean_anomaly

planet1 = ttvfast.models.Planet(
    mass=0.0001303490487,                         # M_sun
    period=19.23891,              # days
    eccentricity=0.0609,
    inclination=88.982,         # degrees
    longnode=0.,           # degrees
    argument=90.,            # degrees
    mean_anomaly=0.,       # degrees
)

planet2 = ttvfast.models.Planet(
    mass=8.98076785e-5,
    period=38.9853,
    eccentricity=0.06691,
    inclination=89.188,
    longnode=0.,
    argument=90.,
    mean_anomaly=0.,
)

planets = [planet1, planet2]
Time = 0                                    # days
dt = 1e-3                                       # days
Total = 1000                                    # days

results = ttvfast.ttvfast(planets, stellar_mass, Time, dt, Total)#, input_flag=0)

#Planet1

idx_planet1 = np.array(results['positions'][0]) == 0
planet1_epochs = np.array(results['positions'][1])[idx_planet1]
planet1_times = np.array(results['positions'][2])[idx_planet1]

#make epoch matrix
A = np.vstack([np.ones(planet1_epochs.size), planet1_epochs]).T

#clean -2 from arrays
cond = planet1_times == -2. 
planet1_epochs = planet1_epochs[~cond]
planet1_times = planet1_times[~cond] #np.where(~cond, planet1_times, 1)# 
A = np.copy(A[~cond])
c, m = np.linalg.lstsq(A, planet1_times, rcond=-1)[0] #leastsq to get best t0 and p

#Planet2

idx_planet2 = np.array(results['positions'][0]) == 1
planet2_epochs = np.array(results['positions'][1])[idx_planet2]
planet2_times = np.array(results['positions'][2])[idx_planet2]

A2 = np.vstack([np.ones(planet2_epochs.size), planet2_epochs]).T

#clean same as before
cond2 = planet2_times == -2.
planet2_epochs = planet2_epochs[~cond2]
planet2_times = planet2_times[~cond2]
A2 = np.copy(A2[~cond2])

c2, m2 = np.linalg.lstsq(A2, planet2_times, rcond=-1)[0] 

fig, ax = plt.subplots(1,1)
ax.plot(A.T[1], (planet1_times-m*A.T[1]-c)*(60.*24.), '.', label='Kepler-9b')
ax.plot(A2.T[1], (planet2_times-m2*A2.T[1]-c2)*(60.*24.), '.', label='Kepler-9c')
ax.set_xlabel("Transit number")
ax.set_ylabel("TTV [min]")
ax.legend()

Questions (0) Why do I constantly get -2, is there a way to simulate the TTVs where I do not get it. From the c version README "BAD_TRANSIT = -1 /* If the bisection method is passed a window that does not bracket the root uniquely, or if the bisection method does not converge within MAX_ITER iteractions, the code returns -1, as discussed in the paper", However I get -2 not -1. Not sure if I should worry about that. In the code above I remove these -2 making sure I do not lose track of the transit times.

(1) If I run the code for different Total = 250, 500, 1000, days I get a different TTV amplitude and shape. I was expecting that a larger Total time would only give me more TTV points, therefore for a zoom at x = [0,12] I should have the same TTV amplitudes more or less. Why not? The figures from top to bottom are for total 250, 500, 1000. I did some tests changing dt to shorter integration steps, and also initial times to negative values.

(2) Usually we do not have the following parameters longnode, argument, mean_anomaly. If I want to reproduce observed TTVs, how exactly may I set these parameters by hand. In the example above, I believe I set both planets to start integration at periapse passage. Moreover, as in question (1) If I change the initial orbital configuration, the TTV shape changes drastically (depending on the values) do you know why? I cannot get a TTV simulation close to the amplitudes from the paper as seen in the image.

I really appreciate any comment/help, thank you!

[simonrw]

Hi @douglasalvesastro12, I'm afraid I don't think I can be of much help, sorry. The python wrapper around ttvfast is pretty faithful to the original (from what I remember), and I'm not too familiar with how it works. I would suggest trying to get the C version working on your input, and maybe taking it up with the original author in this repo: https://github.com/kdeck/TTVFast.

[douglasalvesastro12]

@mindriot101 it's alright, I will contact the author there to get some help. It just got to my mind that perhaps doing a chi2(some_params) I may be able to get some agreement at least. If I get something interesting I may post it in case someone experience similar issues. Anyway, thanks @mindriot101