A Simulation of Transit Timing Variation

Ziyue Zeng

Guangdong Experimental High School, Guangzhou, China.

DOI: 10.4236/jamp.2019.78127                

Abstract

Exoplanet transit timing variation is a method to find exoplanets. To understand this method better, I wrote a computer program in python to simulate the transit of exoplanets. I use my program to simulate the transit timing variation observed in the Kepler-19 system. I make a simple simulation of Kepler-19 system, and this simulation shows that the variation in transit timing due to other planets is very obvious for Kepler-19b, the transiting planet, which means the transit timing variation method is very useful for finding exoplanet in Kepler-19 system. The whole paper is an illustration for that. The simulation I make is relatively simple but it still corresponds to the law of TTV, and because of its simplicity, it can help more people understand.

Keywords

Exoplanet, Transit Timing Variation, Transit Timing, Kepler-19 System

1. Introduction

This paper presents a program that relies primarily on transit photometry and transit timing in the search for exoplanets to simulate the transit timing variation observed in the Kepler-19 system. Transit timing variation is a phenomenon that the period of the planet transiting the star varies due to the gravitational effects of other planets orbiting the star. And we can use this phenomenon to discover other planets which don’t transit base on the transit photometry method, which means that the light curve of a star can indicate the existence of an exoplanet that does not transit. There are also other methods to discover exoplanet indirectly, and transit timing variation and other methods can check each other to make sure an exoplanet really exists. So far, the transit timing variation method has only discovered two exoplanets, Kepler-19c [1] and Kepler-9d. Before 2012, the radial-velocity method was the most productive technique for finding exoplanets. After 2012, the transit method becomes the most productive because of the use of the Kepler spacecraft. But due to their limitations, many exoplanets will be ignored. For example, if the orbital plane of an exoplanet is perpendicular to the line connecting the star surrounded by this exoplanet to the Earth, then we will not be able to observe the change in radial velocity and cannot directly know the existence of this exoplanet through transit, which means both radial-velocity method and transit photometry method don’t work in this case. However, transit timing variation can solve this problem and that’s why I am so interested in it. TTV is very useful as a supplement and I want to know about its feasibility and want to know whether it can do more than being a supplement. What’s more, I hope that more people can know more about it. I wrote a computer program to simulate the transit timing variation of Kepler-19b to better understand the system. The program uses the initial velocity to update force and acceleration, and then uses force and acceleration to update positioning. Current positioning then updates the velocity, and the process repeats. I calculated the initial velocity using Kepler’s Second Law and Kepler’s Third Law and I calculated the difference between periods to make a period variation graph to show the transit timing variation of Kepler-19b due to Kepler-19c and Kepler-19d. I use geometry to calculate the area of the star that is covered by the transiting planet, and this can be used to draw a light curve of Kepler-19.

2. Simulations

2.1. Updating the Position with an Initial Velocity

Following are the equations I use to update the position of the first planet.

$v_{\text{final}}=v_{\text{initial}}+a\times t$ (1)

$v_{\text{final}}=v_{\text{initial}}+\frac{F\times t}{m}$ (2)

$x_{\text{final}}=x_{\text{initial}}+v_{\text{final}}\times t$ (3)

By expressing and calculating velocity, force, and position in terms of x, y, and z coordinates, we can plot them in a three-dimensional coordinate system. In order to use Equation (2) above to update the velocity, we need to calculate the forces between the planet and the star using Newton’s Gravitational law. G is the gravitational constant, M is the mass of the star, and m is the mass of the planet.

$\stackrel{\to }{F}=\frac{G\times M\times m}{r^2}\hat{r}$ (4)

To use Equation (4) and to express and calculate force in terms of x, y, and z direction, we need to calculate the Euclidean distances in x, y, and z direction.

The force in x direction is given by

$F_x=F\times \frac{r_x}{r}\times \hat{r}$

Similarly, the force in y and z direction can be calculated using the method above.

As a result, we can calculate the velocity in terms of x, y, and z direction and then update the position in terms of x, y, and z direction.

We can use the new position to calculate the new forces between celestial bodies and then use the new forces to update the velocity and position in x, y, and z direction again.

In this process, knowing the mass of the celestial bodies, the only variation that I need is the initial velocity and the initial Euclidean distance between the two celestial bodies.

All the methods mentioned above allow us to simulate the motion of many celestial bodies. To get a graph about the motion of many celestial bodies, we can record every position that the celestial bodies pass through.

I randomly set some properties of celestial bodies to check whether my simulation can work (see Table 1).

In order to make the planets orbit around the star, we need to know something about the initial velocity:

If

$v_{\text{initial}}=\sqrt{\frac{GM}{r}\times \hat{r}}\mathrm{,}$

planets will have a circular orbit.

If

$\sqrt{\frac{GM}{r}\times \hat{r}}<v_{\text{initial}}<\sqrt{\frac{2GM}{r}\times \hat{r}},$

planets will have an elliptical orbit.

If

$v_{\text{initial}}=\sqrt{\frac{2GM}{r}\times \hat{r}},$

planets will have a hyperbolic orbit and it will escape from the gravitational field of the star.

The formulas for explanation are as follows:

The orbital velocity:

$F_c=F_g$

$\frac{m{v}_1^2}{r}=\frac{GMm}{r^2}$

$v_1=\sqrt{\frac{GM}{r}}$

The escape velocity:

$KE=U$

$\frac{1}{2}m{v}_2^2=\frac{GMm}{r}$

$v_2=\sqrt{\frac{2GM}{r}}$

And then I make a graph to check.

Table 1. Some random data of celestial bodies. Celestial body 1 is the star. Initial velocities in y and z directions are 0.

Figure 1. The application of the process describe above.

2.2. Using Maths to Calculate the Initial Velocity

In reality, most of the orbits of celestial bodies are ellipses instead of circle. And the central celestial bodies will at one focus of the ellipse. To calculate the initial velocity of the celestial bodies, we need to use some formulas of ellipse.

The semi-major axis is a, the semi-minor axis is b, and the focal length is c.

Eccentricity:

$e=\frac{c}{a}$

semi-minor axis:

$b^2=a^2-c^2$

Area of ellipse:

$\text{area}=\text{π}ab$

For a line between the star and the planet, the line will sweep out an area as the planet moves around the star for a small angle $\text{d}\theta$ the area will be:

$\text{d}A=\frac{1}{2}\times {\left(a-c\right)}^2\times \text{d}\theta$

As a result, the initial velocity can be calculated by combining the formulas above using Kepler’s second law:

$\frac{\text{d}A}{\text{d}t}=\frac{1}{2}\times {\left(a-c\right)}^2\times \frac{\text{d}\theta }{\text{d}t}=\frac{\text{π}ab}{\text{period}}$

$\frac{1}{2}\times \left(a-c\right)\times v_{\text{initial}}=\frac{\text{π}ab}{\text{period}}$

$v_{\text{initial}}=\frac{2\text{π}a\times \sqrt{a^2\times \left(1-e^2\right)}}{\text{period}\times \left(a-e\times a\right)}$

Actually, we can apply all the things mentioned in this section on celestial bodies which have circular orbit. The eccentricity of a circle is 0 and the semi-major axis is the radius of a circle. The formula of calculating the initial velocity of a circular-orbit planet will be:

$v_{\text{initial}}=\frac{2\times \text{π}\times \text{radius}\times \sqrt{{\text{radius}}^2\times \left(1-0^2\right)}}{\text{period}\times \left(\text{radius}-0\times a\right)}$

$v_{\text{initial}}=\frac{2\times \text{π}\times \text{raidus}}{\text{period}}=\frac{\text{perimeter}}{\text{time}}$

So we only need the semi-major axis, eccentricity, and the period to calculate the initial velocity.

2.3. Period Variation

The period of a planet is not constant because of the gravitational effects due to other planets. The period variation is one kind of transit timing variation.

In order to know the differences between periods, I plot periods in a graph in a relatively long time. The way I calculate one period for a particular planet is to record every time when one of the planet passes from negative to positive. And then I calculate the differences between every two record time, which show us the periods of the planets. I randomly choose some data to draw the graph. See Table 2.

2.4. Drawing a Light Curve

To draw a light curve of the star, we need the relationship between the area blocked by the transiting planet and the surface area of the star that we can observe when it is not covered by the planet. To calculate the relationship between these two areas, we need to establish a two-dimensional system and put the center of the star at the origin. Let the distance between the center of the planet and the center of the star be d (see Figure 2) $d={\text{position}}_{\text{star}}-{\text{position}}_{\text{planet}}$ . Because the radius of the orbit of the planet is relatively small compared to the distance between the whole system and the observer, we can let the radius of the planet be $R_p$ and let the radius of the star be $R_s$ . Also because the arc is relatively small compared to the perimeter of the orbit, we can consider the velocity that planet moves toward the star as the linear velocity of the planet. The method of calculating the star’s area that is blocked by the planet is changing with respect to the change in the relative position of the planet and the star.

Figure 2. Period variation of a planet due to the gravitational effects of other planets.

Table 2. Some random properties of celestial bodies. Celestial body 1 is the star.

There are six conditions that need different methods to calculate the blocked area shown in Figure 3.

1) when the planet is not blocking the star: area = 0;

2) when $d\le -\left(R_s+R_p\right)$ and $d\ge -R_s$ ;

3) when $d<-R_s$ and $d\ge d>-\left(R_s-R_p\right)$ ;

4) when $d>-\left(R_s-R_p\right)$ and;

5) when and;

6) when and.

An example of calculating the blocked area is as follow.

In condition 2 shown in Figure 4, the blocked area = area 1 + area 2. We can calculate the sectorial area of the star and the planet first. And then we can calculate area1 by subtract the sectorial area of the star with triangle OAB, and use the same method to calculate area 2.

We can use some similar methods to calculate the areas in other conditions.

Using different methods in these conditions, I get the light curve of the star (see Figure 5).

2.5. Using the Data of Kepler-19

I use the data of Kepler-19 of Wikipedia (shown in Table 3) to run my program assuming that all the orbits of the Kepler-19 system are in the same plane. Following are my results (see Figures 6-8).

Table 3. The data of Kepler-19 system.

Figure 3. A transiting planet.

Figure 4. The way to calculate blocked area in condition 2.

Figure 5. The change in brightness of the star. x axis is time in seconds and y axis is the ratio of the brightness of the star.

Figure 6. The position of the celestial bodies in Kepler-19 system.

Figure 7. The variation of period.

Figure 8. The change in brightness of Kepler-19. x axis is time in seconds and y axis is the ratio of the brightness of the star.

3. Conclusion

I use python to simulate the gravitational effects on celestial bodies due to other celestial bodies, and I get the orbit of the planets orbiting around Kepler-19 and the light curve of Kepler-19. However, giving the limitation of my laptop, the time step I use is not short enough to get accurate simulations. Using the time step of 100 seconds and a total time of 360,000,000 seconds, it takes 1 hour and 20 minutes to run the program. To get a more accurate simulation, I need a better computer to run the program. What’s more, I will try to do a program that can analyze the light curve of a star to get the information about the mass, radius, and many other properties of the planets which transits can be observed by us and of all the planets that are orbiting around the star.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]

Ballard, S., Fabrycky, D., Fressin, F., Charbonneau, D., Desert, J.-M., Torres, G., Marcy, G., Burke, C.J., Isaacson, H., Henze, C., Steffen, J.H., Ciardi, D.R., Howell, S.B., Cochran, W.D., Endl, M., Bryson, S.T., Rowe, J.F., Holman, M.J., Lissauer, J.J., Jenkins, J.M., Still, M., Ford, E.B., Christiansen, J.L., Middour, C.K., Haas, M.R., Li, J., Hall, J.R., McCauliff, S., Batalha, N.M., Koch, D.G. and Borucki, W.J. (2011) The Kepler-19 System: A Transiting 2.2 R⊕ Planet and a Second Planet Detected via Transit Timing Variations. The Astrophysical Journal, 743, 200. https://doi.org/10.1088/0004-637X/743/2/200