Simulations for Terrestrial Planets Formation
In this paper, we investigate the formation of terrestrial planets in the late stage of planetary formation using two-planet model. At that time, the protostar has formed for about 3 Myr and the gas disk has dissipated. In the model, the perturbations from Jupiter and Saturn are considered. We also consider variations of the mass of outer planet, and the initial eccentricities and inclinations of embryos and planetesimals. Our results show that, terrestrial planets are formed in 50 Myr, and the accretion rate is about . In each simulation, 3 - 4 terrestrial planets are formed inside ”Jupiter” with masses of . In the 0.5 - 4 AU, when the eccentricities of planetesimals are excited, planetesimals are able to accrete material from wide radial direction. The plenty of water material of the terrestrial planet in the Habitable Zone may be transferred from the farther places by this mechanism. Accretion could also happen a few times between two major planets only if the outer planet has a moderate mass and the small terrestrial planet could survive at some resonances over time scale of yr. In one of our simulations, com-mensurability of the orbital periods of planets is very common. Moreover, a librating-circulating 3:2 configuration of mean motion resonance is found.
Subject headings:astrophysics-exoplanet-planetary formation-n-body simulation
Since the discovery of the first extrasolar planet around Solar-Type star, the detection of the extrasolar planets develops rapidly. To date, more than planets are found orbiting their center stars beyond our solar system, including 35 multiple planetary systems. Recently, scientists have found evidences of methane and carbon dioxide in the atmosphere of a Hot-Jupiter (HD 189733 b) (http://planetquest.jpl.nasa.gov). Of known extrasolar planets, the minimum mass is generally several Jupiter masses. There are also several terrestrial planets (Super-Earth), but their orbital characteristic is unsuitable for the formation or development of life. Along with the development of survey techniques and incoming high definition space missions, people will definitely discover more and more Earth-like planets in the extrasolar planetary systems. The research of formation and evolution of the terrestrial planet now becomes important topics in astrophysics, astrobiology, astrochemistry and so on.
Planet formation has certain order (Zhou et al., 2005), and Jupiter-like planets at greater distance are formed faster than those near the Sun. It is generally believed that the planet formation may experience the following stages: The grains condensed in the initial stage grow to km-sized planetesimals in the early stage, and then, in the middle stage, Moon-to-Mars sized embryos are formed by accretion of the planetesimals. The size of embryos correlates with the feeding zone of the planetesimals. According to the formula of Hill radius: (where , are the heliocentric distance and the mass of planetesimal), the distant planetesimals have wider feeding zones, so the formed embryos are larger. When the embryos grow to the core of mass ( ), runaway accretion may take place accordingly. With more atmosphere accreted, the embryos contract, growing ever denser and more massive, eventually collapse to form giant Jovian planets (Hu & Xu, 2008; Ida & Lin, 2004). However, at the same period, the inner planetesimals accrete in respective accretion scope, and then the embryos of terrestrial planet (namely kind of the terrestrial planet core matter) are formed. At the end of the third stage, it is around that the protostar has formed for about Myr, the gas disk has dissipated. A few larger bodies with low and are in crowds of planetesimals with certain eccentricities , and inclinations . In the late stage, the terrestrial planetary embryos are excited to high eccentricity orbit by gravitational perturbation. Then, the orbital crossing makes the planets accreting material in the broader radial area. Solid residue is either scattered out of the planetary system or accreted by the massive planet. However, it also has the possibility of being captured at the resonance position of the major planets (Nagasawa & Ida, 2000; Hu & Xu, 2008).
Taking our Solar System as the background, Chambers (2001) made a study of terrestrial planet formation in the late stage by numerical simulations. He set Moon-to-Mars size planetary embryos in the area of AU, include gravitational perturbations from Jupiter and Saturn. He also examined two initial mass distributions: approximately uniform masses, and a bimodal mass distribution. The results show that planets are formed in Myr, and finally survive over Myr timescale. The space distribution and concentration (see Section 4 in Chambers (2001)) of planets formed in the simulations are similar to our solar system. However, the planets produced by the simulations usually have eccentric orbits with higher eccentricities , and inclinations than Venus. Raymond, Quinn & Lunine (2004, 2006) also studied the formation of terrestrial planets. In the simulations, they took into account Jupiter’s gravitational perturbation, wider distribution of material ( AU) and higher resolution. The results confirm a leading hypothesis for the origin of Earth’s water: they may come from the material in the outer area by impacts in the late stage of planet formation. Raymond, Mandell, & Sigurdsson (2006) explored the planet formation under planetary migration of the giant. In the simulations, super Hot Earth form interior to the migrating giant planet, and water-rich, Earth-size terrestrial planet are present in the Habitable Zone ( AU) and can survive over yr timescale.
In our model, Solar System is taken as the background. But several changes are worth noting : 1) we use two-planet (Jupiter and Saturn, see Fig.1) model. 2) in the model, Jupiter and Saturn are supposed to be formed at the beginning of the simulation, with two swarms of planetesimals distributed among AU and AU respectively. 3) The initial eccentricities and inclinations of planetesimals are considered. 4) The variations of the mass of Saturn are examined. 5) The exchange of material in the radial direction is also studied by the parameter of water mass fraction. 6) We perform the simulations over longer timescale ( Myr) in order to check the stability and the dynamical structure evolution of the system. Our results show that the terrestrial planets produced interior to Jupiter have higher mass accretion rate, and share the similar architecture as the Solar System. However, the structure beyond Jupiter correlates with the initial mass of Saturn. Almost each simulation has a water-rich terrestrial planet in the Habitable Zone ( AU).
In Section 2, the initial conditions, algorithm and integration procedure are described in detail. Section 3 presents the main results. We conclude the outcomes in Section 4.
2.1. Initial conditions
Generally, the time scale for formation of Jupiter-like planet is less than Myr (Briceño, 2001). Nevertheless, as we know, the formation scenario of planet embryos is related to their heliocentric distances and the initial mass of the star nebular. If we consider the model of MMSN (minimum mass solar nebular), the upper bound of the time scale for Jupiter-like planet formation corresponds to the time scale for the embryo formation at AU (Kokubo & Ida, 2002), which is just at resonance location of Jupiter. In the region AU, embryos will be cleared off by strong gravitational perturbation arising from Jupiter. There should be some much smaller solid residue among Jupiter and Saturn, even though the ’clearing effect’ may throw out most of the material in this area. That’s why we set embryos only in the region AU and planetesimals in the AU and AU.
We adopt the surface density profile as follows (Raymond, Quinn & Lunine, 2004):
In (1), is the surface density at snowline, the snowline is at AU with .As mentioned earlier, the mass of planetary embryos is proportional to the width of the feeding zone, which is associated with Hill Radius, , so the mass of an embryo increases as
The embryos among AU are spaced by ( varying randomly between 2 and 5) mutual Hill Radii, , which is defined as
where and are the semi-major axes and masses of the embryos respectively. Replacing in (2) with , and substituting (1) in (2), then, we get relations between the mass of embryos and the parameter as
As shown in Fig. 1, the initial planetesimals are spread over AU (excluding Hill Radii around Jupiter); the distribution of them should meet (1). Here, we equally set the masses of planetesimals inside and outside Jupiter respectively as shown in (5). Consequently, the number distribution of the planetesimals is only needed to satisfy . Additionally, we keep the total number of planetesimals and embryos inside Jupiter, and the number of planetesimals outside Jupiter both equal to .
The water mass fraction of the bodies is same as Raymond, Quinn & Lunine (2004), i.e., the planetesimals beyond 2.5 AU have water material by mass, those between AU have water material by mass, and the others have water material by mass. The eccentricities and inclinations vary in () and (), respectively. The mass of Saturn in simulations 1a/1b,2a/2b and 3a/3b are , , respectively. Each simulation is carried out twice with a) considering, b) not considering self-gravitation of planetesimals among Jupiter and Saturn.
In regular coordinate system, the motion equations of an n-body system are (Murray & Dermott, 1999)
where the index denotes the body , and , are the generalized coordinate and momentum of the body , respectively. The Hamiltonian, , is the sum of the kinetic and potential energy for the system. From (6), we know that the rate of any quantity, , can be conveniently expressed in the following form,
If we define an operator (Chambers, 1999), then we can rewrite (7) as . Integral the differential equation over time (), and so we get
where and are the values of corresponding to the time and respectively. If we define ( actually is the internal time step), then expand (8) at zero, we have
The symplectic integrator is to divide into pieces, each piece could be individually solved, and then they approximate the solution of the problem via applying the solutions once a time. For example, we split the Hamiltonian , and hence have the operators . It is easy to obtain from (8). By expanding the exponential (attn. ), ignoring the second- and higher-order small quantities of , then we have
We can get a second-order integrator by applying a small equivalence transformation.
The key point for symplectic algorithm is how to split Hamiltonian into pieces. Considering a dynamical system composed of bodies orbiting a massive central body, we can split the Hamiltonian into the primary and the secondary parts. Chambers (1999) proposed a hybrid symplectic algorithm, in which the Hamiltonian is divided into the following parts:
where is the unperturbed Keplerian motion of the smaller bodies, is the total interaction potential energy of the smaller bodies, and is the kinetic energy of the center body (Note That: has different meaning from above, refers to the numbers of bodies excluding the central body). The term of ’hybrid’ means that, for the convenience of calculation, heliocentric coordinates and barycentric velocities are used while solving (11) (Chambers, 1999; Duncan, Levison & Lee, 1998). From (10), we can infer that all of the higher terms depend on both and . If (), and therefore, the second-order integrator is correct to only when the different parts of Hamiltonian meet the conditions . However, when close encounter occurs between two bodies, the distance between them approaches zero, hence the can not be satisfied. Chambers (1999) introduced a changeover function to translate part of associated with the close encounter to , and then integrate it using Bulirsch-Stoer method (Stoer & Bulirsch, 1980). The modified are present as
tends to zero when is small, while tending to one when is large (Chambers, 1999).
We use the hybrid symplectic integrator (Chambers, 1999) in MERCURY package to integrate all the simulations. We take into account that collision and coalescence will occur, when the minimum distance between any of the two objects is equal to or less than the sum of their physical radii. While they were separated by not more than 3 Hill radii, we consider close encounters will take place. When the distance from the central star is more than AU, these bodies are removed, because they are so far from the central star that they play an insignificant role of the interaction. In addition, we adopt days as the length of time step, which is a twentieth period of the innermost body at AU. The simulations are carried out over Myr time scale. At the end of the intergration, the changes of energy and angular momenta are and respectively. The simulations are performed on a workstation composed of CPUs with GHz, and each costs roughly days.
All of simulations exhibit some classical processes on planet formation. Firstly, we will analyse simulation 2a/2b to discuss the physical processes which can apply to every simulation. Next, we will make a statistical analysis in order to find out that how the planet formation may rely on different physical factors.
3.1. Simulation 2a/2b
At the end of the calculation, terrestrial planets are formed in 2a/2b. Table 1 shows the properties of the terrestrial planets from simulations 2a and 2b. We label the planets as b, c, d, and so on, according to the heliocentric distance (hereinafter). The masses of the terrestrial planets range from several Mars masses to several Earth masses. All of them are water-rich, except the planet e in simulation 2b. Some parameters of a certain planet are comparable with the terrestrial planets in solar system. For example, the orbital eccentricity of planet b in simulation 2b is , which is very close to that of Earth.
Fig. 2 is a snapshot of simulation 2a. At Myr, it is clear that the planetesimals are excited at the ( AU), ( AU) and ( AU) resonance locations with Jupiter, and this is similar to the Kirkwood gaps of the asteroidal belt in solar system. For about Myr, planetesimals and embryos are deeply intermixed, most of the bodies have large eccentricities. Collisions and accretions emerge among planetesimals and embryos. This process continues until about Myr, the planetary embryos are mostly formed, and then dynamical evolution is start. The formation time scale in our work is in accordance with that of (Ida & Lin, 2004). Finally, inside Jupiter, terrestrial planets are formed with masses of . However, at the outer region, planetesimals are continuously scattered out of the system at Myr. For about Myr, there are no survivals except at some resonances with the giant planet. As shown in Fig. 2, there is a small body at the resonance with Jupiter. Due to the planetesimals’ scattering, Jupiter (Saturn) migrates inward (outward) AU ( AU) toward the sun respectively. Such kind of migration agrees with the work of Fernandez et al. (Fernandez & Ip, 1984) Hence, the mean motion resonance is destroyed, then the ratios of periods between Jupiter and Saturn degenerate to . Therefore, the ratio of periods for Jupiter, small body and Saturn is approximate to .
Fig. 3 is a snapshot of simulation 2b. In comparison with Fig. 2, it is apparent that planetesimals are excited more quickly at the ( AU), ( AU) and ( AU) resonance location with Jupiter. The several characteristic time scales are the same as simulation 2a for the bodies within Jupiter. 4 planets are formed in simulation 2b, the changes of position of Jupiter and Sat-urn are about the same as simulation 2a. We note simulations 2a and 2b have the same initial conditions, the only difference between them is whether we consider the self-gravitation among the outer planetesimals. The results of simulation 2b are shown to be a consequence of being expected but not surprising. There is a little stack planetesimals survival over Myr among AU, located in the area of ( AU) and ( AU) resonances with Jupiter.
Planets in 55 cnc planetary system have similar spatial distribution to the solar system (Fischer et al., 2008), From Table 1 and Fig. 3, it is not difficult to see that 4 terrestrial planets formed in simulation 2b move on the nearly-circular orbit. Mars is ever regarded as a survivor of an original planetary embryo, according to its unique chemical and isotopic characteristics. As a matter of fact, the planet e in simulation 2b is a survivor of the initial planetary embryos. Planet e in simulation 2b does not accrete anything over Myr integration time. Comparing the semi-major axis of Mars with that of planet e, we notice that the planetesimals of simulation 2b are located in the asteroidal belt. Furthermore, the ratio of periods of the planet e and Jupiter is nearly . The total mass of the main-belt in solar system is about , being (Hu & Xu, 2008) of the initial solid material. If the assumed planet e in simulation 2b would break into thousands of fragments, they may undergo re-accretion or ejection by the perturbation of Jupiter over secular evolution. If it is similar to the same ratio of mass of the belt in solar system, then the leftovers of the solid materials almost bear a total mass of . In this sense, an asteroidal belt is very likely to form in the system quite similar to that of our solar system.
Fig. 4a is the mass curve of terrestrial planets for simulation 2b. We can find that the accrete velocity is not uniform. At Myr, planets reach half mass of their final mass, and then, the accretion velocity slows down, because the planetesimals are only a quarter left. Until about Myr, the terrestrial planets are formed. The accretion rate and mass concentration (the mass rate of the largest terrestrial planet and the total formed objects) are and , respectively. The corresponding parameters in simulation 2a are and respectively. However, in the area outside Jupiter, initial material is scattered out of the system.
Planet embryos are formed from feeding zones where the planetesimals are located. A feeding zone has unique chemical and isotopic characteristics. It is helpful to study the trace of the planetesimals to understand the composition of terrestrial planet, vice versa, for example, if we can investigate the origin and formation process by revealing the chemical or isotopic characteristics of the moon. In Fig. 4b is shown the trace of survivals. We note that all the materials accreted by terrestrial planets come from the inner swarm of planetesimals or embryos. Here Jupiter is like a wall, which separates the inner and outer planetesimals from exchanging materials. Once again, Fig. 4b verifies that Jupiter may protect the inner terrestrial planets from colliding with the outer bodies (Wetherill, 1990). We set a water mass fraction on each body, it is easy to work out how much water-material of the finally terrestrial planet bears. Take planet c in simulation 2b for example, the water material is approximately , about 8 times than Earth. Fig. 4b can also verify that terrestrial planet accrete material in broad radial direction.
3.2. Statistical analysis
The production efficiency of the terrestrial planet in our model is high, and the accretion rate inside Jupiter is in the simulations. terrestrial planets formed in Myr. 5 of 6 simulations have a terrestrial planet in the Habitable Zone ( AU) (see Fig. 5). The planetary systems are formed to have nearly circular orbit and coplanarity, similar to the solar system (see Table 2). We suppose that the above characteristics are correlated with the initial small eccentricities and inclinations. Such adoption could generate more close encounters or collisions in the early several Myr, which may increase viscosity of the system and then make the orbits more nested on circular orbit on the orbital plane. The concentration in Table 2 means the ratio of maximum terrestrial planet formed in the simulation and the total terrestrial planets mass. It represents different capability on accretion, and is not associated with self-gravitation. The average value of this parameter is similar to the solar system. Considering the self-gravitation of planetesimals among Jupiter and Saturn, the system has a better viscosity, so that the planetesimals will be excited slower. The consideration of self-gravitation may not change the formation time scale of terrestrial planets, but will affect the initial accretion speed and the eventual accretion rate. In Table 2, the simulations 1b, 2b, 3b have a bit higher accretion rate. When the self-gravitation is not considered, the planetesimals may be excited quickly. The accretion has a faster speed at the early several Myr, so this can promote the accretion rate.
From Fig. 4b, we have to be aware of that Saturn accretes a few planetesimals, this is uncommon in simulations 1a, 1b, 3a, 3b. And Fig. 5 shows the finally structure of the simulations, it is clear that different Saturn mass will affect the outer structure of the system beyond Jupiter. In simulation 3a (3b), the Saturn mass is . Now it is large enough to clear the area among Jupiter and Saturn. In simulation 1a (1b), Saturn’s mass is , more or less equal to the embryos’ mass. In this case, it is too small to clear off any planetesimal amongst the region of Jupiter and Saturn. Therefore, we draw the conclusion that only the Saturn’s mass is close to be , then accretion may happen. Saturn and Jupiter in our solar system may form in same stage. If there exist embryos of outside Saturn, the giant Jupiter-mass planets may form.
The scattering of planetesimals could cause the migration of planets, for example, Jupiter migrated from AU to AU while Saturn traveled from AU to AU. There are plenty of ratios of semi-major axis of the survival planets nearly . In this case, the planet is easily to be captured on mean motion resonance (Lee & Peale, 2002; Lee, 2004; Zhou et al., 2005; Fischer et al., 2008). There are still some ratios of periods between the survival bodies close to a simple ratio of integers (see Fig. 5). In the very long process of dynamical evolution after planetary formation, the planets also have the possibility of been captured onto resonant orbit. For example, the orbits show the mean motion resonance be-tween Jupiter and outer small body in simulation 1b, and the resonance between Jupiter and the planet d in simulation 3a and so on. Many researchers have studied the resonance and stability of the planetary systems (Ji et al., 2003; Zhou & Sun, 2003; Zhou et al., 2004). As shown in Fig. 6, in simulation 2a, two planets are on crossing orbits. When a close encounter occurs, mean motion resonance is formed, with resonance angle (where are the mean longitude and the longitudes of periapse, the footnotes 1, 2 means the inner and outer planets respectively.) librating around , while circulating, and shows large oscillations. Such ’librating-circulating’ configuration is similar to the configuration of resonance in HD 73526 planetary system. There have been several hypotheses about its origin (Tinney et al., 2006; Sándor & Kley, 2006; Sándor, Kley & Klagyivik, 2007). However, it still needs further study in the future.
We simulate the terrestrial planets formation by using two-planet model. In the simulation, the variations of the mass of outer planet, the initial eccentricities and inclinations of embryos and planetesimals are also considered. The results show that, during the terrestrial planets formation, planets can accrete material from different regions inside Jupiter. Among AU, the accretion rate of terrestrial planet is , i.e., about initial mass is removed during the progress. The planetesimals will improve the efficiency of accretion rate for certain initial eccentricities and inclinations, and this also makes the newly-born terrestrial planets have lower orbital eccentricities. It is maybe a common phenomenon in the planet formation that the water-rich terrestrial planet is formed in the Habitable Zone. The structure, which is similar to that of solar system, may explain the results of disintegration of a terrestrial planet. Most of the planetesimals among Jupiter and Saturn are scattered out of the planetary systems, and this migration caused by scattering (Fernandez & Ip, 1984) or long-term orbital evolution can make planets capture at some mean motion resonance location. Accretion could also happen a few times between two planets if the outer planet has a moderate mass, and the small terrestrial planet could survive at some resonances over yr time scale. Structurally, Saturn has little effect on the architecture inside Jupiter, owing to its protection. However, obviously, a different Saturn mass could play a vital role of the structure outer Jupiter. Jupiter and Saturn in the solar system may form over the same period.
In our simulations, neither terrestrial planets are formed within AU, nor planetesimals or embryos are left. However, a lot of exoplanets with orbital semi-major less than AU are observed, and several Super-Earths are discovered. It is usually believed that they were formed far from the center star and then migrated into current location (Raymond, Mandell, & Sigurdsson, 2006). We do not consider the migration in the simulations, which is caused by the interaction between the giant planets or planetesimals in the gaseous disk (Ou et al., 2007). So the simulation in this work can be applied to the case of the dissipation of gas disk, in the late stage of planet formation. In the future study, we will consider the giant planets under inward migration, and in such circumstances short-period terrestrial planet could be produced. To date, terrestrial planets are not detected in the observations, due to the reasons of the selection effect of detection methods and the low resolution precision. Both Doppler velocities and transit method are sensitive to the objects moving in smaller orbits. The current research of extrasolar terrestrial planets has greatly contributed to the origin and evolution of our own solar system. Kepler has been launched successfully on March 6, 2009, whose main scientific objective is detecting the Earth-like terrestrial planets. Along with high accuracy incoming space projects, it is predictable that more and more extrasolar planetary systems with similar structure to the solar system will be discovered.
- affiliation: Graduate School of Chinese Academy of Sciences, Beijing 100049, China
- affiliation: Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, China
- affiliation: Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, China
- affiliation: National Astronomical Observatory, Chinese Academy of Sciences,Beijing 100012, China
- Briceño, C. et al., 2001, Science, 291, 93
- Chambers, J. E. 1999, MNRAS, 304, 793
- Chambers, J. E. 2001, Icarus, 152, 205
- Duncan, M., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067
- Fernandez, J. A., & Ip, W. H. 1984, Icarus, 58, 109
- Fischer, D., et al. 2008, ApJ, 675, 790
- Hu, Z. W., & Xu, W. B. 2008, Planet Science, Beijing: Science Press (in Chinese)
- Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388
- Ji, J. H., et al.2003, ApJ, 585, L139
- Ji, J. H., et al. 2003, ApJ, 591, L57
- Ji, J. H., et al. 2003, Chin. Astron. Astrophysics, 27, 127
- Kokubo, E., & Ida, S. 2002, ApJ, 581, 666
- Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596
- Lee, M. H. 2004, ApJ, 611, 517
- Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics, New York: Cambridge University Press
- Nagasawa, M., & Ida, S. 2000, AJ, 120, 3311
- Ou, S. L., et al. 2007, ApJ, 667, 1220
- Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icarus, 168, 1
- Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265
- Raymond, S. N., Mandell, A. M., & Sigurdsson, S. 2006, Science, 313, 1413
- Sándor, Z., & Kley, W. 2006, A&A, 451, L31
- Sándor, Z., Kley, W., & Klagyivik, P. 2007, å, 472, 981 28
- Stoer, J., & Bulirsch, R. 1980, Introduction to Numerical Analysis, New York: Springer-Verlag
- Tinney, C. G., et al. 2006, ApJ, 2006, 647, 594
- Wetherill, G. W. 1990, Ann. Rev. Earth Planet Sci., 18, 205
- Zhou, L. Y., et al. 2004, MNRAS, 350, 1495
- Zhou, J. L., & Sun, Y. S. 2003, ApJ, 598, 1290
- Zhou, J. L., et al. 2005, ApJ, 631, L85