rodluger / planetplanet

A general photodynamical code for exoplanet light curves
https://rodluger.github.io/planetplanet
GNU General Public License v3.0
45 stars 1 forks source link

What does a transit time really mean? #56

Closed rodluger closed 7 years ago

rodluger commented 7 years ago

Currently I am drawing from the following distributions (mu, sigma) for the transit times of each of the TRAPPIST-1 planets:

transits = [(7322.51736, 0.00010), 
            (7282.80728, 0.00019), 
            (7670.14165, 0.00035), 
            (7660.37859, 0.00038),
            (7671.39767, 0.00023), 
            (7665.34937, 0.00021), 
            (7662.55284, 0.00037)]

I took these from Gillon et al. (2017) and Luger et al. (2017). Notice, however, that the transit times for b and c are defined for transits that occur several hundred days before the those of the outer planets. This is a problem when running an N-body code because if I start my integration around day 7322 (close to t0 for planet b), there will be sufficient time for planet-planet interactions to shift the times of transit of the other planets away from the value t0 we specified. In other words, even though I'm setting t0 = 7662.55284 for planet h, if I start the integration several hundred days before that, the planet will likely not transit at that time.

This issue is rooted in how I use t0 to specify the orbits: I use it to compute the time of pericenter passage:

fi = (3 * np.pi / 2.) - w
tperi0 = t0 + (per * np.sqrt(1. - ecc * ecc) / (2. * np.pi) * (ecc * np.sin(fi) /
              (1. + ecc * np.cos(fi)) - 2. / np.sqrt(1. - ecc * ecc) *
              np.arctan2(np.sqrt(1. - ecc * ecc) * np.tan(fi / 2.), 1. + ecc)))

which in turn is used to compute the mean anomaly:

M = 2. * PI / per * modulus(time[0] - tperi0, per);

The mean anomaly is used to get the eccentric anomaly and then the true anomaly at t = 0, which is finally passed to REBOUND to initialize the orbits. The problem is in the step where I compute the mean anomaly: I'm assuming the planet is on a perfectly periodic orbit, which allows me to take the modulus.

Currently, this doesn't matter much, since we're only doing PPO statistics. But if we eventually want to use this to predict and model PPOs, we need our transit time priors to be physically meaningful. The only way I see to do this correctly is to run some sort of MCMC chain every time we draw from the prior to iteratively adjust tperi0 for each planet so that they are actually seen to transit at their respective t0s.

@ericagol: What's the best approach here? I bet this is a common problem when doing TTV analyses. Perhaps linear/analytic TTV theory could come to the rescue here?

ericagol commented 7 years ago

The best way to handle this would be to draw from the TTV model posterior.

However, the TTV posterior is continuing to improve due to additional data, and so I don't think this is necessary for this paper.

I would recommend choosing the transit times closest to the starting time of the integration (these times are given online with a link to the file in the relevant figure of the Nature paper). That way the TTVs won't have changed much for the individual planets between the starting point of the N-body & the time of first transit.

Sent from my iPad

On Aug 25, 2017, at 2:04 AM, Rodrigo Luger notifications@github.com wrote:

Currently I am drawing from the following distributions (mu, sigma) for the transit times of each of the TRAPPIST-1 planets:

transits = [(7322.51736, 0.00010), (7282.80728, 0.00019), (7670.14165, 0.00035), (7660.37859, 0.00038), (7671.39767, 0.00023), (7665.34937, 0.00021), (7662.55284, 0.00037)] I took these from Gillon et al. (2017) and Luger et al. (2017). Notice, however, that the transit times for b and c are defined for transits that occur several hundred days before the those of the outer planets. This is a problem when running an N-body code because if I start my integration around day 7322 (close to t0 for planet b), there will be sufficient time for planet-planet interactions to shift the times of transit of the other planets away from the value t0 we specified. In other words, even though I'm setting t0 = 7662.55284 for planet h, if I start the integration several hundred days before that, the planet will likely not transit at that time.

This issue is rooted in how I use t0 to specify the orbits: I use it to compute the time of pericenter passage:

fi = (3 np.pi / 2.) - w tperi0 = t0 + (per np.sqrt(1. - ecc ecc) / (2. np.pi) (ecc np.sin(fi) / (1. + ecc np.cos(fi)) - 2. / np.sqrt(1. - ecc ecc) np.arctan2(np.sqrt(1. - ecc ecc) * np.tan(fi / 2.), 1. + ecc))) which in turn is used to compute the mean anomaly:

M = 2. PI / per modulus(time[0] - tperi0, per); The mean anomaly is used to get the eccentric anomaly and then the true anomaly at t = 0, which is finally passed to REBOUND to initialize the orbits. The problem is in the step where I compute the mean anomaly: I'm assuming the planet is on a perfectly periodic orbit, which allows me to take the modulus.

Currently, this doesn't matter much, since we're only doing PPO statistics. But if we eventually want to use this to predict and model PPOs, we need our transit time priors to be physically meaningful. The only way I see to do this correctly is to run some sort of MCMC chain every time we draw from the prior to iteratively adjust tperi0 for each planet so that they are actually seen to transit at their respective t0s.

@ericagol: What's the best approach here? I bet this is a common problem when doing TTV analyses. Perhaps linear/analytic TTV theory could come to the rescue here?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

rodluger commented 7 years ago

Awesome! I don't see a file with the posteriors, but they did have a table with all the Spitzer transit times. I am now using the transit times closest to BJD - 2,450,000 = 7670 (the last transits they report), which is on October 8, 2016. I'm adding comments throughout the code that our ephemerides are only accurate close to that date, and that we will update them once we have new observations.

On Sat, Aug 26, 2017 at 1:18 AM Eric Agol notifications@github.com wrote:

The best way to handle this would be to draw from the TTV model posterior.

However, the TTV posterior is continuing to improve due to additional data, and so I don't think this is necessary for this paper.

I would recommend choosing the transit times closest to the starting time of the integration (these times are given online with a link to the file in the relevant figure of the Nature paper). That way the TTVs won't have changed much for the individual planets between the starting point of the N-body & the time of first transit.

Sent from my iPad

On Aug 25, 2017, at 2:04 AM, Rodrigo Luger notifications@github.com wrote:

Currently I am drawing from the following distributions (mu, sigma) for the transit times of each of the TRAPPIST-1 planets:

transits = [(7322.51736, 0.00010), (7282.80728, 0.00019), (7670.14165, 0.00035), (7660.37859, 0.00038), (7671.39767, 0.00023), (7665.34937, 0.00021), (7662.55284, 0.00037)] I took these from Gillon et al. (2017) and Luger et al. (2017). Notice, however, that the transit times for b and c are defined for transits that occur several hundred days before the those of the outer planets. This is a problem when running an N-body code because if I start my integration around day 7322 (close to t0 for planet b), there will be sufficient time for planet-planet interactions to shift the times of transit of the other planets away from the value t0 we specified. In other words, even though I'm setting t0 = 7662.55284 for planet h, if I start the integration several hundred days before that, the planet will likely not transit at that time.

This issue is rooted in how I use t0 to specify the orbits: I use it to compute the time of pericenter passage:

fi = (3 np.pi / 2.) - w tperi0 = t0 + (per np.sqrt(1. - ecc ecc) / (2. np.pi) (ecc np.sin(fi) / (1. + ecc np.cos(fi)) - 2. / np.sqrt(1. - ecc ecc) np.arctan2(np.sqrt(1. - ecc ecc) * np.tan(fi / 2.), 1. + ecc))) which in turn is used to compute the mean anomaly:

M = 2. PI / per modulus(time[0] - tperi0, per); The mean anomaly is used to get the eccentric anomaly and then the true anomaly at t = 0, which is finally passed to REBOUND to initialize the orbits. The problem is in the step where I compute the mean anomaly: I'm assuming the planet is on a perfectly periodic orbit, which allows me to take the modulus.

Currently, this doesn't matter much, since we're only doing PPO statistics. But if we eventually want to use this to predict and model PPOs, we need our transit time priors to be physically meaningful. The only way I see to do this correctly is to run some sort of MCMC chain every time we draw from the prior to iteratively adjust tperi0 for each planet so that they are actually seen to transit at their respective t0s.

@ericagol: What's the best approach here? I bet this is a common problem when doing TTV analyses. Perhaps linear/analytic TTV theory could come to the rescue here?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/rodluger/planetplanet/issues/56#issuecomment-325101718, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK_sTuwDNJEXIXQEICvEnO7_Cxajdks5sb9TogaJpZM4PCE-J .

rodluger commented 7 years ago

I'm closing this for now. We'll revisit it in the future.