A&A 463, 359-367 (2007)
DOI: 10.1051/0004-6361:20066582

Co-orbital terrestrial planets in exoplanetary systems: a formation scenario

C. Beaugé1 - Zs. Sándor2 - B. Érdi2 - Á. Süli2


1 - Observatorio Astronómico, Universidad Nacional de Córdoba, Laprida 854, (X5000BGR) Córdoba, Argentina
2 - Department of Astronomy, Loránd Eötvös University, Pázmány Péter sétány. 1/A, 1117 Budapest, Hungary

Received 17 October 2006 / Accepted 11 November 2006

Abstract
Aims. We study the formation of a hypothetical terrestrial-type body in the equilateral Lagrange points of a giant extrasolar planet. Starting from a swarm of planetesimals in stable tadpole orbits, we simulate its dynamical and collisional evolution under a wide range of different initial conditions and masses for both the Trojan population and its planetary companion. We also analyze the effects of gas drag from the interaction of the planetesimals with the nebular disk.
Methods. The formation process is simulated with an N-body code that considers full gravitational interactions between the planetesimals and the giant planet. Gas interaction is modeled with Stokes and Epstein drags, where the drag coefficients are chosen following the results of full hydrodynamic simulations performed with the 2D public hydro-code FARGO.
Results. In both gas-free and gas-rich scenarios, we have been able to obtain a single final terrestrial-type body in a stable tadpole orbit around one of the triangular Lagrange points of the system. However, due to gravitational instabilities within the swarm, the accretional process is not very efficient and the mass of the final planet never seems to exceed $\sim$0.6 Earth masses, even when the total mass of the swarm is five times this value. Finally, we also included an orbital decay of the giant planet due to a type II migration. Although the accretional process shows evidence of a lower efficiency, a small terrestrial planet is still able to form, and follows the giant planet towards the habitable zone of the hosting star.

Key words: celestial mechanics - planets and satellites: formation - methods: N-body simulations

1 Introduction

After the discovery of the first extrasolar planet around 51 Pegasi by Mayor & Queloz (1995), the most important milestone of extrasolar planetary research would certainly be the discovery of the first Earth-sized planet. Several studies have appeared in the recent years concerning the dynamics and stability of these still hypothetical Earth-like planets. The ultimate aim of these works and the extrasolar planetary research is to find analogy between our Solar System and the other extrasolar planetary systems (EPS), and therefore to understand the formation and evolution of the planetary systems in a more global way.

There are several studies dealing with a rather specific class of hypothetical Earth-like planets which are co-orbital companions of a gas giant. The possible existence of such a Trojan-like planet in the habitable zone of an EPS has been mentioned first by Menou & Tabachnik (2003). The first detailed study of dynamical habitability of particular systems has been performed by Dvorak et al. (2004) and Érdi & Sándor (2005). There are other investigations as well not excluding the case of a more massive Trojan planet, and not restrictive to the habitable zones (see Gozdziewski & Konacki 2006; Laughlin & Chambers 2002; Nauenberg 2002). These studies are also based on a Solar System analogy. It is well known that around the stable Lagrangian points of the Sun-Jupiter system there exists a large population of small celestial bodies known as the Trojan asteroids. The dynamics of Trojans has always been a very fascinating subject of celestial mechanics, several authors studied the dynamical particularity of their motion (e.g. Érdi 1997, for a review), the size of the stability region (Levison et al. 1997; Efthymiopoulos & Sándor 2005; Giorgilli & Skokos 1997), and their formation (see Morbidelli et al. 2005; Marzari & Scholl 1998).

Works related to the extrasolar Trojans (or exotrojans), however, assume that the three-body system is already in existence. Thus, studies of exotrojans have so far concentrated on the dynamical stability as a function of the orbital parameters of the known giant planet, and/or the mass of the Earth-type body. In this paper, however, we will focus our attention on a possible formation mechanism for such a hypothetical body. In particular, we wish to address two points: (i) under what circumstances a planetesimal swarm orbiting around one of the equilateral Lagrange points will accrete into a terrestrial mass body? (ii) If the giant planet and the associated co-orbital body formed far from the star, what is the effect of the inward migration process on the stability of the Earth-like Trojan planet?

2 Physical scenario for the formation of the co-orbital system

According to current theories of planetary formation, terrestrial planets formed through accretional collisions of rocky planetesimals. We thus begin our study assuming the existence of a swarm of planetesimals trapped in the stability region around a stable Lagrangian point (L4 or L5) of the star-giant planet system. This population may have been captured by any of several mechanisms, such as gas drag in the co-orbital region, mass growth of the planetary companion, etc. The giant planet itself is supposed to be located far from the star having a small orbital eccentricity.

In the course of their orbital evolution, the planetesimals will collide and, if their relative velocities are sufficiently small, larger bodies may be formed. Having larger mass, these embryos may enhance the relative velocities of the rest of the swarm, favoring further collisions, ending in the possible formation of a single rocky planet. The key point of this scenario is that the relative velocities should be small enough avoiding the escape of the planetesimals after the collision events. A mechanism which can reduce the relative velocities of the planetesimals could be the collision itself. However, if the co-orbital region is still embedded in the gaseous protoplanetary disk, gas drag may help the accretion process.

Since we assume an inward migration of the giant planet, be it simultaneous or posterior to the collisional evolution of the swarm, we will also incorporate a low-density gaseous disk enveloping the complete system. In recent years, the effects of the so-called type II migration on three-body systems has been investigated by several authors using gravitational three-body calculations (e.g. Beaugé et al. 2006; Lee & Peale 2002) and full hydrodynamical simulations (e.g. Masset 2001; Kley 2003). In the case of mean-motion resonances of type pn1 - qn2=0 (where ni are the orbital frequencies and $p \ne q$), numerical integrations indicate that the commensurability is maintained during migration. Although the same not been proved for co-orbital bodies, we can hope that the Trojan planetesimals will remain locked in the 1:1 resonance, and an exotrojan planet formed far from the star may reach the habitable region in a stable orbit.

The presence of a nebular gas does not necessarily imply a gas-rich formation scenario, and the environment in the vicinity of the equilateral Lagrange points depends on the mass of the giant planet. If it is smaller than a critical value $M_{\rm c}$ which depends on the disk's parameters (e.g. Ward 1997; Crida et al. 2006), both the giant planet and the co-orbital regions are embedded in the gas, and its dynamical effects on the planetesimals must be taken into consideration. However, for planetary masses larger $M_{\rm c}$, the giant planet opens a gap in the disk (e.g. Lin & Papaloizou 1980, 1993; Kley 1999), and the gas density in the co-orbital regions reaches values close to zero. In such a case, the dynamical evolution of a Trojan population will be practically gas-free, even if a significant nebular disk is still present in the exoplanetary system.

In our simulations, we will analyze both the gas-free and gas-rich scenarios. We will also consider cases where the inward migration of the giant planet occurs after, or simultaneous with, the formation process of the exotrojan body.

3 Evolution through collisional accretion without gas drag

A gas-free scenario can be the result of an efficient gap opening mechanism, if the formation timescale of a Trojan-like planet is much larger than the gap opening timescale. However, the same scenario can also describe the case where the exotrojan population was trapped in co-orbital motion after the evaporation of the gas material by the central star. In our simulations, the mass of the hosting star M* was taken equal to one Solar mass, and the giant planet M1 was considered to be a Jupiter clone, having the same mass and present osculating orbital elements. However, some runs were also performed with other planetary masses and eccentricities. A typical initial planetesimal population consisted of 500 equal-mass bodies with a total mass of  $1.5~M_{\oplus}$. Their initial conditions were chosen as follows:

\begin{displaymath}\begin{array}{lll}
a \in [5.1,5.3] \; {\rm AU} \;\;\; & e \in...
...80] \; {\rm deg} & \Omega \in [0,360] \;
{\rm deg}
\end{array}\end{displaymath} (1)

where a is the semimajor axis, e the eccentricity, I the inclination, M the mean anomaly, $\varpi$ the longitude of pericenter and $\Omega$ the longitude of the ascending node. These initial conditions place the bodies deep inside the libration region of the L4 equilateral Lagrange point of the system. Additional runs were also performed in L5.

The swarm was allowed to evolve due to the gravitational interactions with the star, giant planet and the other members of the population. We used a Bulirsch-Stoer integrator with variable step size. Whenever the mutual distance of two planetesimals mi and mj was smaller than five times their mutual physical radii, they were accreted, creating a new body with a mass equal to mi+mj and its new orbit was given by the center of mass approximation. Throughout the simulations we assumed a material density of 2 g/cm3. Due to the use of an inflated radius for each planetesimal, the accretion rate will be accelerated by a factor somewhere between f or f2 (Kokubo & Ida 1996), where f is the ratio between the adopted and real physical radius (in our case f=5). Thus, it is important to keep in mind that the formation timescales obtained in all our simulations cannot be taken on face value, and must be scaled accordingly.

We performed a series of runs according to this scenario, adopting different values for the total number of bodies N0 (from 100 to 1000) and total mass $M_{\rm T}$ (from 1 to  $3~M_{\oplus}$). Computational resources defined the limit on N0, since we consider full mutual interactions on all planetesimals. However, it must be recalled that we are working inside a relatively small area of the phase space, restricted to a portion of the tadpole region of the L4 point. Even with this consideration, we still expect some numerical dispersion in our results. Its magnitude was checked performing a series of simulations for each value of N0 and $M_{\rm T}$, with different random values for the orbital elements of the swarm, and comparing the results. We found that the final mass of the formed Trojan planet did not depend significantly neither on the initial mass, nor on the number of the whole population of the planetesimals. We also performed numerical simulations for smaller mass giant planet ( $30{-}100~M_{\oplus}$) but the average mass of the Trojan planet was $\sim$ $0.3~M_{\oplus}$ never exceeding  $0.6~M_{\oplus}$. (We note that in the case of smaller mass giant planet the width of the stability region around L4/L5 scales roughly as $\sim$ (m1/(m0+m1))1/2, thus we took into account this fact by decreasing the radial extent of the initial population.)

The results of two typical runs are shown in Figs. 1 and 2. The first corresponds to an initial population of N0=1000 bodies and a total mass of $M_{\rm T}=3~M_{\oplus}$, while the second has a smaller swarm (N0=300) and half the total mass. In each case, the top plot shows the temporal variation of three parameters: at a given time, N indicates the number of surviving bodies of the swarm, $n_{\rm col}$ is the cumulative number of accretional collisions and  $n_{\rm esc}$ gives the cumulative number of ejections of planetesimals from the system. Ejections occurred whenever one of the bodies escaped from the tadpole region of the Lagrange point and suffered a close encounter with the giant planet. The middle plots show the mass  $m_{\rm max}$ of the largest body in the swarm, as well as the mean value of the mass $\langle m \rangle$. Finally, the bottom frames give the evolution of the mean eccentricity $\langle e \rangle$ and mean inclination (in radians) $\langle I \rangle$ of the Trojan population.


  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{6582fig1.eps}
\end{figure} Figure 1: Typical simulation of the accretion process of a Trojan-like planet in a gas-free scenario. Size of the initial population was N0=1000 and total mass was $M_{\rm T}=3~M_{\oplus}$. Top plot shows the number of surviving planetesimals N as function of time, as well as the cumulative number of collisions $n_{\rm col}$ and escapes $n_{\rm esc}$. The middle graph shows the evolution of the average mass of the swarm $\langle m \rangle$ and it largest value  $m_{\rm max}$, while the bottom plot gives the average eccentricity $\langle e \rangle$ (black) and inclination $\langle I \rangle$ (grey). Recall that we are using inflated physical radii for the planetesimals, so the time axis must be scaled accordingly.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{6582fig2.eps}
\end{figure} Figure 2: Same as previous figure, but with N0=300 and $M_{\rm T}=1.5~M_{\oplus}$. Note that the time axis is now in linear scale.
Open with DEXTER

Qualitatively, both simulations show the same global behavior. Accretion starts soon after the start of the runs and, as the final outcome of the dynamical and collisional evolution of the system, yields a single massive body in a stable tadpole configuration. A first stage of the simulation is characterized by all the members of the swarm well restricted to the tadpole region. Collisions rapidly give rise to a marked mass spectrum, where runaway accretion causes the appearance of a few (usually one or two) planetary embryos but with most of the rest of the population still in the original mass range. In a second stage, however, several escapes were noted. In Fig. 1 this occurs around $t \sim 10^3$ years, while in Fig. 2 it happens later, at $t \sim 10^4$ years (recall that the times must be scaled according to the inflated radii). Although the first cases of ejections occur for very small planetesimals, as the system evolved we also find ejections of more massive bodies. Around $t \sim 5 \times 10^4$ years, this gravitational instability seemed to dominate the evolution of the system. There were practically no more accretions, and the additional decrease in the number of bodies is mainly caused by ejections. This evaporation of the remnant population takes longer in Fig. 1 than in Fig. 2, but the behavior seems very similar. Thus, the size of the biggest surviving Trojan-like object appears to be defined at some mid point in the formation process, and its final mass is only equal to approximately  $0.4~M_{\oplus}$.

The bottom frames of both figures show that the eccentricity of this final body oscillates around 0.05, a value very similar to that of the planet. Even though this was the type of result we hoped for, it appears that the accretional process is not very efficient, and at least two thirds of the initial population was lost. In fact, even though the initial mass in both runs differed by a factor of two, the final mass of the exotrojan planet was very similar. A comparison with additional runs showed that the first escapes always seemed to occur when the mass of the biggest embryo reached a value of $\sim$ $0.1{-}0.15~M_{\oplus}$, although the time associated with the first ejections could vary. Moreover, even though the total initial mass varied between 1 and $3~M_{\oplus}$, the mass of the final planet was at most $\sim$ $0.6~M_{\oplus}$.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{6582fig3.eps}
\end{figure} Figure 3: Osculating eccentricity and semimajor axis of the Trojan swarm in four different time frames during the accretional process. The area of the circles are proportional to the planetesimal's mass.
Open with DEXTER

Figure 3 shows the eccentricity and semimajor axis distribution of the Trojan swarm at four different times during the evolution of the simulation shown in Fig. 2. The size of the circles is proportional to the planetesimal's mass. Larger bodies show a tendency to remain near the center of the librational region, while smaller planetesimals have a larger dispersion.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{6582fig4.eps}
\end{figure} Figure 4: The size of the stability region as a function of the mass of a Trojan planetesimal placed to L4 calculated by the RLI chaos indicator. Light colors indicate stable, while dark colors indicate unstable orbits. It can be seen that a $0.15~M_{\oplus}$ Trojan planetesimal can seriously destabilize the region around L4.
Open with DEXTER

This dynamical instability of a Trojan body undergoing the additional perturbation of a massive planetesimal inside the Lagrange stability region was previously unknown. In order to quantify this effect, we performed additional numerical investigations by using the relative Lyapunov indicator (RLI) chaos detection method (Sándor et al. 2004). We placed a massive embryo exactly in the L4 point, whose mass ($M_{\rm p}$) was varied between $0{-}0.15~M_{\oplus}$ with $\Delta M_{\rm p}=0.002~M_{\oplus}$. For each value of $M_{\rm p}$ we determined the width of the stability region placing massless Trojan bodies equidistantly from 4.9-5.5 AU with $\Delta a=0.004$ AU on the line connecting the star to the Lagrangian point L4. The giant planet M1 was once again taken to be a Jupiter clone. Results are shown in Fig. 4, where the horizontal axis shows the semi-major axis of the massless Trojans, while the vertical axis corresponds to the mass $M_{\rm p}$ of the massive Trojan embryo. We note that a $M_{\rm p} > 0.15~M_{\oplus}$ Trojan planetesimal can destabilize almost the whole Trojan population, while stable orbits within the co-orbital region are only found for smaller masses. We should note that the above stability study represents the smallest perturbation for the Trojan swarm since the massive Trojan is placed exactly in L4. If this was deviated from L4 it would librate thus representing a stronger perturbation to the Trojan swarm.

Figure 5 shows similar stability plots, only this time we show the maximum eccentricity attained by each initial condition over a time interval T. White points correspond to initial conditions were the change in eccentricity was small (stable), and dark points correspond to increases in eccentricity of the order of unity (unstable orbits). The total integration time was T=103 orbital periods (top graph) and T=104 orbital periods (bottom frame). From the last graph it is clear that, for long time spans, even small values of $M_{\rm p}$ are able to generate instabilities in all of the tadpole region.

We repeated these calculations for smaller values of M1, under the assumption that the accretional process may have began during the mass growth phase of the giant planet. In all cases, however, the same type of behavior was noted, and the instability actually begin sooner for even smaller values of $M_{\rm p}$.


  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{6582fig5.eps}\par\vspace*{2mm}
\includegraphics[width=7.8cm,clip]{6582fig6.eps}
\end{figure} Figure 5: Same as previous figure, now showing the logarithm of the maximum eccentricity attained by each initial condition during the complete integration timespan T. Top: T = 103 orbital periods. Bottom: T = 104 orbital periods. Note that the region in black corresponds to a change in the eccentricity of the order of unity; consequently the orbit is unstable over the timespan T.
Open with DEXTER

This gravitationally instability as a stopping mechanism for the formation process shows an interesting contrast with the formation of terrestrial planets from planetesimals in (usual) heliocentric orbits. In both cases, the appearance of planetary embryos induce chaotic and unstable motion of the small bodies in their vicinity. However, in the case of heliocentric orbits, this orbital instability actually helps the accretional process since it causes an excitation of the eccentricities and helps to avoid the isolation of the feeding zones of each embryo. Planetesimals on unstable heliocentric orbits have no where to go, except continue to cross the orbit of the embryos until they are finally accreted. The case of a swarm in the vicinity of the co-orbital region is different. The stable tadpole zone has a limited extent, and unstable orbits do not remain in its vicinity, but are subjected to close encounters with the giant planet, resulting in the ejection from the system.

Thus, in the Trojan scenario, the planetary formation process itself is responsible for its own ulterior inefficiency. It is therefore doubtful whether a final planet larger than a few tenth of  $M_{\oplus}$ is possible, since the step from $\sim$ $0.1~M_{\oplus}$ onwards is done in an ever increasing hostile environment. Summarizing our simulations in a gas-free scenario we can conclude that the formation of a Trojan-like planet around the stable Lagrangian points is possible through pure collisional accretion mechanism, although the final mass always seems to be smaller than one Earth-mass. In the next section we shall investigate the formation mechanism helped by the presence of gas material.

4 Gas-rich collisional accretion

The origin of the Jovian Trojan asteroids in our own Solar System is not yet well established. Classical theories place the capture of these bodies during the formation stage of the outer planets. According to this school, planetesimals in heliocentric orbits in the vicinity of the growing Jupiter may have been captured into stable tadpole orbits either by drag effects of the gaseous nebula (e.g. Shoemaker et al. 1989) or the rising of the potential peaks at the triangular points due to Jupiter's mass growth (e.g. Marzari & Scholl 1998). A different version was suggested by Peale (1993), according to which physical collisions among small planetesimals, either in temporary horseshoe or nearby heliocentric orbits, may have lead to the formation of larger bodies in stable librations. Note that this scenario is similar to the one discussed previously in this work, the main difference being the final outcome of the collisional evolution.

However, recent simulations by Morbidelli et al. (2005) have considered the possibility that the present Trojan population may have been captured during the orbital migration suffered by the outer planets during the last stages of planetary formation. It is well known (Michtchenko et al. 2001; Gomes 1998) that a passage of Jupiter and Saturn through the 2/1 mean-motion resonance destabilizes all Trojan orbits of the Jovian system. According to Morbidelli et al. (2005), the same mechanism that could have caused the ejection of any primordial Trojans, may also have opened the road to temporary captures of new bodies. As both giant planets left the resonant configuration, the Trojan region would regain its orbital stability, and any new transitory members would become permanent.

Such a chaotic capture mechanism requires the existence of at least two giant planets in nearby orbits and a past temporary passage through a low order mean-motion resonance. None of the extrasolar planetary systems that are the ultimate goal of this work belong to this category. Most of them only contain one known giant planet or, if any companion is also observed, it is located far from any important resonant configuration. Of course, this does not rule out the existence of undetected giant planets in these systems, but such an analysis is beyond the scope of the present work. Consequently, it appears at least possible that any hypothetical exotrojans in these planetary systems may have been captured at a time when the nebular gas was still present and, perhaps, the giant planet was still in its formation stage.

With such a scenario in mind, the main objective of this section is to analyze whether a growing giant planet and/or dynamical interactions between the Trojans and the gas may affect the accretional evolution of the swarm. Both planet growth of the giant planet and gas drag are believed to cause a damping effect on the tadpole orbits (e.g. Fleming & Hamilton 2000), and could thus help to stabilize the otherwise gravitationally unstable population and increase the efficiency of the accretional process. Laughlin & Chambers (2002) also analyzed the gas dynamics of a nebular disk in the vicinity of L5, and observed a certain accumulation of gaseous material in its vicinity. This led the authors to postulate that a Trojan-type planet may be formed in such a configuration. However, it is important to keep in mind that the presence of gas around the Lagrange points is only temporary, and is emptied depending on the viscosity as soon as the gap opening by the giant planet is complete (Kley, private communication). Moreover, as we will show presently, the main effects of the gas on a solid planetary core is gas drag.

Our first task in this regard is to model the gas-Trojan interactions. Past works on this topic have used gas drag models, such as Stokes or Epstein approximations (e.g. Marzari & Scholl 1998), or the application of the adiabatic invariant theory (Sicardy & Dubois 2003; Fleming & Hamilton 2000). Although these approaches may be appropriate for our own Solar System, they may be too approximate in extrasolar planetary systems where the disk-planet interactions were much more intense up to a point of leading to a large scale planetary migration.

In order to derive and test an adequate model for the effects of gas dynamics in the orbital evolution of tadpole orbits, we first performed a series of full hydrodynamical simulations using the FARGO 2D public code developed by F.S. Masset, and that can be found in http://www.maths.qmul.ac.uk/~masset/fargo/. This code solves the Navier-Stokes equations for a Keplerian disk subject to the gravity of the central object and that of embedded planets. The disk is isothermal and not self-gravitating. The algorithm behind the code can be found in Masset (2000).

  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{6582fig7.eps}\par\vspace*{2mm}
\par\includegraphics[width=8cm,clip]{6582fig8.eps}
\end{figure} Figure 6: Hydrodynamical simulations of a gaseous disk with two embedded planets, performed with FARGO. Top: gas density plots. Darker (lighter) regions correspond to smaller (larger) surface density values. Left plot corresponds to a giant planet with mass $M_1=30~M_{\oplus}$, while the right graph was obtained with $M_1=300~M_{\oplus}$. In both frames the location of the giant planet is shown with a white circle in the X axis, while the position of the Trojan body can be seen in the vicinity of the L4 point. Bottom: surface density of the gaseous disk, as function of the distance from the star, for an azimuthal angle crossing the L4 point. The three curves correspond to simulations with masses of the giant planet equal to 30, 100 and  $300~M_{\oplus}$.
Open with DEXTER

We placed a giant planet with the same initial conditions as in the first section (except that the semimajor axis was taken as unity) plus a small Trojan-like body with mass equal to  $10^{-3}~M_{\oplus}$. Both were placed in elliptic orbits with eccentricity 0.05. While the Trojan body was always allowed to feel the gravitational effects of the gas, in this series of runs the giant planet was not affected. This was done in order to avoid secular changes in the planets' orbital elements which could introduce additional dynamical effects on the orbital motion of the small planetesimal. Finally, the gaseous disk was assumed to have an $\nu$-viscosity equal to 10-5, constant surface density $\sigma_0 = 6.4 \times 10^{-4}$ and extended from 0.4 to 2.5 AU.

Figure 6 (top) presents snapshots of two simulations, showing the levels of surface density of the gas in grey scale. Lighter (darker) shades correspond to higher (lower) density values. The position of the giant planet is shown with a white circle in the X axis, while the small Trojan-like body is located in the vicinity of the L4 equilateral point. When the giant planet's mass is smaller than a critical value $M_{\rm c}$, a significant amount of gas is still present in the Trojan region. However, when the mass is large enough to generate a gap in the disk, the Lagrange point appears empty. The lower plot of Fig. 6 shows the variation of the gas density as a function of the distance to the star, and calculated along a radial direction which crosses the L4 point. For a small planet with $M_1 \sim 30~M_{\oplus}$, the gas density in the Trojan region is practically the same as in the other regions of the disk. Thus, we expect that the Trojan-disk interactions should be significant. Larger masses of the giant planet open a gap in the disk, which includes the Trojan tadpole region. The density of the gas in the orbital region of the Trojan is thus very diminished. For giant planets in the Jovian mass range the gap-opening appears complete and the orbital evolution of the Trojan should follow similar roads as in a gas-free scenario.

It is known (e.g. Ward 1997; Crida et al. 2006) that the critical mass for gap opening is a function of the disk parameters, such as the $\alpha$-viscosity and the relative disk scale height H/r, and can also be dependent on the physics underlying the hydrodynamical model. Thus, although the three "phases'' shown on the lower frame of Fig. 6 should occur at some time during the mass growth of the giant planet, they do not necessarily have to be tied to these mass values. However, it must be recalled that our current aim is not to model the planet-disk interactions, but simply to qualitatively analyze what is the dynamical effect on a Trojan body in each case. Figure 7 shows the evolution of the orbital eccentricity (left plots) and librational amplitude (right) of the test Trojan for each different value of M1. The first case corresponds to a gas-rich environment, and leads to a damping of the eccentricity towards a value $e \sim 0.02$. The second case corresponds to $M_1 = 100~M_{\oplus}$. Some gas is still present in the tadpole region, although the gas density is not homogeneous and some time-dependent spiral structures can be observed in the area. The eccentricity still suffers a damping effect and, unexpectedly, the asymptotic value is lower than for higher gas density. Finally, the lower frames in Fig. 7 correspond to a complete gap opening and practically no gas in the co-orbital region. Consequently, the eccentricity shows little secular effect, although we do note an increase in the amplitudes of oscillation. In all cases, however, the libration amplitude of the critical argument of the 1:1 resonance shows no appreciable damping. On the contrary, for a small planetary mass there actually seems to be a small but noticeable increase in the amplitude.


  \begin{figure}
\par\includegraphics[width=8.2cm,clip]{6582fig9.eps}
\end{figure} Figure 7: Black: results of three hydrodynamical simulations of a system consisting of a massive planet of mass m1 and a small Trojan-like body near the L4 Lagrange point. Left frames show the evolution of the Trojan's orbital eccentricity, while the right frames correspond to the amplitude of libration $\lambda - \lambda _1$. Top graphs were obtained with a massive planet of mass $M_1=30~M_{\oplus}$, middle graphs with $M_1 = 100~M_{\oplus}$, and bottom plots with $M_1=300~M_{\oplus}$. Grey: temporal variation of the eccentricity and libration amplitude obtained from a 3-body model where the gas effects on the Trojan body are modeled with Stokes gas drag with  C = 10-5.
Open with DEXTER

Knowing the dynamical effects of the Trojan-disk interactions for different initial conditions and several values of the planetary mass, we can attempt to model them using a simple gas drag model that can be incorporated into our N-body code. We tested two well known expressions: Stokes and Epstein (or square-velocity) drag laws. The acceleration felt by a solid spherical body in each case is given by:

                $\displaystyle {\rm Stokes}$ : $\displaystyle \;\; \ddot{\vec{r}} = -{6 \pi \mu s \over m_{\rm p}} {\vec{u}} = -C {\vec{u}}$ (2)
$\displaystyle {\rm Epstein}$ : $\displaystyle \;\; \ddot{\vec{r}} = -C_D {\pi s \over 2 m_{\rm p}} \rho_{\rm g} \vert u\vert \vec{u}$  

where C is the Stokes-drag coefficient, CD is a constant coefficient (usually of the order of unity), $m_{\rm p}$ is the mass of the perturbed body, s its physical radius, $\mu$ is the molecular viscosity of the gas and $\rho_{\rm g}$ its density. Notice that Epstein's law depends explicitly and linearly with the gas density, while Stokes' model is independent of this parameter. In both cases,

\begin{displaymath}{\vec{u}} = \dot{\vec{r}} - \eta (\vec{\omega} \times \vec{r})
\end{displaymath} (3)

is the relative velocity between the solid body and the gas. In this equation, $\vec{r}$ is the astrocentric position vector of the perturbed body and  $\vec{\omega}$ is the circular angular velocity at the point, defined by the two-body problem with the star. The constant $\eta$ is the ratio between the gas velocity and the Keplerian velocity and is usually taken as $\eta = 0.995$ (see Beaugé & Ferraz-Mello 1993). Stokes drag is believed to be an adequate model for aerodynamical forces when gas density is very low (Adachi et al. 1976), such as inside a well defined gap in the disk. For large gas densities, a square-velocity (or Epstein drag) is usually preferred.

We performed several tests with both drag laws added to a 3-body system, where both the massive planet M1 and the Trojan body have the same mass parameters and initial conditions as those used in the hydro-simulations. For the $M_1=30~M_{\oplus}$ case, numerical simulations showed that both gas laws are able to reproduce the damping of the eccentricity for adequate values of the drag coefficients. However, for a Trojan in an L4 libration, the Epstein drag caused a swift and large increase in the librational amplitude $\sigma = \lambda - \lambda_1$, such that the body escaped the tadpole region in timescales of the order of a few hundred years. Although Stokes drag also caused an increase in $\Delta \sigma$, the magnitude was barely appreciable for values of $C \sim 10^{-5}$ 1/day, consistent with an equilibrium eccentricity of $e \sim 0.19$. A comparison of both the effects with the upper frames of Fig. 7 shows that Stokes drag (grey curves) seems to be more compatible with the hydrodynamical simulations.

For the intermediate case with $M_1 = 100~M_{\oplus}$, once again Epstein drag caused an increase in the librational amplitude which is not observed in the hydro simulations. Stokes drag, with the same value of C, did not cause a significant eccentricity damping, since the larger planetary mass was able to counteract the dissipative effect and maintain the orbital eccentricity at a value similar to that of the planet (see the grey curves in the middle frames of Fig. 7). The results of FARGO, however, show that the eccentricity should in fact decrease significantly, and a small but appreciable damping in $\Delta \sigma$ can also be noted. While we were not able to reproduce both these effects with an Epstein drag, a similar behavior was possible with Stokes drag, but adopting a value of $C \sim 10^{-4}$ 1/day, about an order of magnitude larger than for the smaller planetary mass. Finally, for $M_1=300~M_{\oplus}$ and practically no gas in the tadpole region, both gas drag laws were adequate for very small values of the coefficients. In particular, the grey curves in the lower frames of Fig. 7 were obtained with Stokes drag and $C \sim 10^{-6}$ 1/day.

In conclusion, compared with the Epstein drag, we find that a Stokes drag model seems to yield dynamical effects which are more compatible with the hydrodynamical simulations. However, both models have serious difficulties in reproducing the FARGO results for intermediate planetary masses, where the dynamics of the gas inside the tadpole region can be complex. For this reason, we only performed numerical simulations of the formation process of our exotrojan planet with constant planetary mass M1. This approximation is equivalent to assuming that the growth time of the giant planet is much larger than the accretional timescale of the Trojan swarm. The main advantage of this approach was that it permitted us to bypass the variation of the drag magnitude with M1. Thus, for each of our runs, and each value of the planetary mass, we adopted a value of C compatible with a FARGO run.


  \begin{figure}
\par\includegraphics[width=8.2cm,clip]{6582fg10.eps}
\end{figure} Figure 8: Typical simulation of the accretion process of a Trojan-like planetary body in a gas-rich scenario. The size of the initial population was N0=300, the total mass was $M_{\rm T}=1.5~M_{\oplus}$ and the initial drag coefficient was set to C=10-5.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.4cm,clip]{6582fg11.eps}
\end{figure} Figure 9: Top: histogram of the mass of the final Trojan body obtained from two different sets of simulations in a gas-rich scenario. The continuous lines correspond to 50 runs with C=10-5, and the broken lines to 50 simulations with a drag five times greater. The mass of the giant planet was set to $M_1=300~M_{\oplus}$, and the initial mass of the swarm was $M_{\rm T}=1~M_{\oplus}$. Bottom: distribution of the orbital eccentricity of the final bodies.
Open with DEXTER

A typical result is shown in Fig. 8. Although the timescale of the formation process is significantly shorter than for the gas-free scenario, there do not seem to be any other significant differences. Gas drag does not seem to be able to counteract the dynamical instability of the tadpole region once a planetary embryo reaches $0.1{-}0.15~M_{\oplus}$, and once again the efficiency of the accretional process is low. The final planetary mass is of the order of $\sim$ $0.33~M_{\oplus}$. The top plot of Fig. 9 shows two histograms of the final planetary mass obtained from two series of 50 simulations. The difference between these simulations was that the randomly chosen initial conditions for the Trojan planetesimals were varied. In all cases the total mass of the system was  $1~M_{\oplus}$, and the mass of the giant planet was set to  $300~M_{\oplus}$. Continuous lines correspond to simulations with C=10-5, and broken lines to simulations with five times this value. Note that there is an appreciable dispersion in the final Trojan mass, but we did not find any case where it surpassed $0.6~M_{\oplus}$. The lower frame of Fig. 9 shows the dispersion in the final orbital eccentricity. In all cases the value of e was very close to the planetary value.

In conclusion, it does not appear that dynamical effects stemming from a gas drag cause significant differences in the accretional behavior of the system. Although we considered different values of the planetary mass (from 30 to  $300~M_{\oplus}$), time-dependent drag coefficients and mass growth for the giant planet, the main results of all our simulations were qualitatively similar. Very small masses of the giant planet usually lead to smaller Trojan bodies, since the stability of the tadpole region was also reduced, and gravitational instabilities caused by the growing embryos set in sooner. For example, our runs with $M_1 \le 30~M_{\oplus}$ did not yield Trojan bodies larger than $\sim$ $0.3~M_{\oplus}$. However, for $M_1 \ge 50~M_{\oplus}$ we did find final masses of the same order as those obtained for fully grown planetary companions.

5 Planetary migration

So far, we have studied the formation scenario of the exotrojan planet considering that both planets are located far from the central star. We will now consider what are the effects of an orbital migration of the giant planet on the accretional process of the Trojan swarm and to the orbital stability of the final terrestrial planet. We imagine that the semimajor axis of the giant planet decays due to a type II migration, in a timescale between 105-107 years (Ward 1997). Thus, the nebular disk is still present in the system, but not in the region of the Trojan population. Either the exotrojan was already formed at this time, or accretion continued in a gas-free environment.

Our first concern was whether the gas dynamics would allow the Trojan population to accompany the giant planet during its migration towards the inner planetary system. We then performed a new series of hydrodynamical simulations with FARGO. Once again, these included a giant planet (present day Jupiter mass and initial orbit) plus an exotrojan in a stable tadpole orbit with amplitude of approximately 10 degrees. Both planets were allowed to feel the gas. Results showed that even in the case of an incomplete gap formation and a strong gas interaction, the exotrojan migrated together with the planet, and the tadpole configuration was maintained (see Fig. 10). This allows us to, once again, model this evolution introducing an ad hoc migration effect into our N-body code (see Kley et al. 2005; Beaugé et al. 2006, for similar modelling).


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{6582fg12.eps}
\end{figure} Figure 10: Decay in semimajor axis of the giant planet (black) and a Trojan body (grey) obtained with a FARGO hydro-simulation of a type II migration. Planetary masses were $M_1=300~M_{\oplus}$ and $M_{\rm p}=0.001~M_{\oplus}$. Note that the semimajor axis of the Trojan body oscillates about the giant planet, and follows the orbital decay.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{6582fg13.eps}
\end{figure} Figure 11: N-body simulation of the accretional evolution of a Trojan swarm coinciding with an inward type II migration of the giant planet. After $t \sim 10^4$ years, only a single exotrojan remains, and follows the giant planet towards a stable tadpole orbit in the habitable region of the system.
Open with DEXTER

We performed several simulations, using different e-folding times $\tau_{\rm a}$ for the semimajor axis decay. The eccentricity of the giant planet was assumed to remain unchanged during the migration. Figure 11 shows typical results of an accretional process in a Trojan swarm occurring simultaneously with an inward migration of the giant planet. The initial planetesimal swarm consisted of N0=300 bodies with total mass equal to one Earth mass, and were placed in a tadpole configuration around L4. Initial conditions were chosen following the same restrictions as for the other simulations in this paper. The giant planet was started with an initial semimajor axis a1=5.2 AU, and migrated to 1 AU in a timescale of almost 105 years. Upon reaching this distance, migration was switched off. Throughout the orbital evolution, the eccentricity and inclination of M1 remained unchanged, and its value can be seen as a horizontal line in the two lower frames of Fig. 11. The orbital elements of the accreting Trojan population are also plotted in each of the graphs. From the top plot we note that the complete swarm follows closely the orbital decay of the giant planet. However, the decrease in the heliocentric distance implies a decrease in the stability region of the tadpole region, and escapes are more frequent than in "static'' simulations. The eccentricities of the planetesimals can reach high values (around 0.3) before being ejected. Nevertheless, the final outcome is very similar to our previous runs. A single exotrojan survives with $M_{\rm p} \sim 0.4~M_{\oplus}$ exhibiting a very stable libration around L4. Note the variation of this planet's eccentricity and inclination in the lower frames for  $t \sim 10^5$ years.

Notice that we have only considered the effects of a planetary migration due to disk-planet interactions, and we have not considered an orbital decay caused by the gravitational interactions of the giant planet with a remnant planetesimal disk. This second type of migration is believed to have played an important role in the giant planets of our own Solar System (see Fernandez & Ip 1984; Hahn & Malhotra 1999), although it seems less important in the case of exoplanets. Nevertheless, it is worthwhile recalling that Beaugé et al. (2002) showed that interactions with a massive planetesimal disk could lead to orbital instabilities in primordial planetary satellites. Considering the sensitivity of the tadpole region to external perturbations, it is not impossible to imagine that similar fly-bys could also affect the stability of a migrating exotrojan body.

6 Conclusions

The basic aim of this work has been to present a simple study of the formation, dynamical evolution, and stability of a terrestrial-type planet in the equilateral Lagrangian points (i.e. the L4 and L5 points) of a giant extrasolar planet. The planet is initially placed far from the star with mass and orbital characteristics similar to Jupiter, and then allowed to migrate inwards towards the habitable region of the stellar system. For Solar-type stars, this region is roughly between 0.7 and 1.3 AU. In the literature it is possible to find several works devoted to analyzing the orbital stability of terrestrial-type exotrojans, in the case where the giant planet is already in the habitable zone. However, there has been no specific attempt to study the formation process of such a planet out of an initial planetesimal swarm, nor analyze the stability of such a planet during the planetary migration process.

We have found that the accretional process of such a system is not efficient enough to form an Earth-massed planet. Since the tadpole region has a stable region of limited extent, any sufficiently large perturbation can introduce significant orbital instabilities into the planetesimal swarm and seriously affect the accretional process. This perturbation can in fact come from the planetary embryos themselves who are part of the accreting population. Thus, we have found that the formation process of a terrestrial planet in L4 or L5 appears mass-limited. Although we have performed a number of simulation with varied initial conditions and different planetary masses, we have been unable to obtain an exotrojan planet with mass larger than  $0.6~M_{\oplus}$. Introduction of gas dynamics, including gas drag, does not seem to modify this picture significantly, and all our simulations gave final masses of the same order of magnitude. It is worthwhile to mention that the dynamical effects of the gas in the tadpole region is complex. Simple drag laws (such as Stokes or Epstein) are not always adequate, especially when the gap opening is in progress. However, we found that numerical simulations with Stokes drag are able to reproduce the the main aspects of the hydro-simulations, although some fine-tuning of the drag coefficient is necessary. Simulations with Epstein drag were not so successful.

Our results, however, support the formation of a Mars-sized, or even bigger planet having mass $\sim$ $0.3{-}0.6~M_{\oplus}$. Thus the discovery of less massive Earth-like Trojan planets by future space missions might be expected. On the other hand, we should notethat the formation of much larger masses seem unlikely within our formation scenario. Such EPSs containing giant planets in the 1:1 resonance (as suggested by Gozdziewski & Konacki 2006, in the cases of HD 73526 and HD 82943) seem difficult to achieve, unless other formation scenarios are assumed, such as gas-instability collapse or capture of already formed planets.

Finally, planetary migration of the giant planet due to disk-planets interactions does not seem to hinder the formation process in the Trojan population, nor does it lead to orbital instabilities in the tadpole region. Thus, if such a terrestrial-type body could have been formed near L4 of a giant planet prior to its migration towards the habitable zone of an extrasolar planetary system, it is possible that it may constitute a planet that could satisfy some conditions of habitability.

Acknowledgements
This work has been supported jointly by the Secretaria de Ciencia y Tecnologia de la Republica Argentina (SECYT) and by the Hungarian National Office of Research and Technology (NKTH) under the grant ARG-4/2004. Zs.S. and B.É. also acknowledge the supports of the Hungarian Scientific Research Fund (OTKA) grant numbers D048424 and T043739.

References

 

Copyright ESO 2007