Stability and Formation of the Resonant System HD 73526
aa r X i v : . [ a s t r o - ph ] J un Astronomy&Astrophysicsmanuscript no. sandoretal c (cid:13)
ESO 2018June 1, 2018
Stability and Formation of the Resonant System HD 73526
Zs. S´andor, , W. Kley, and P. Klagyivik Department of Astronomy, E¨otv¨os Lor´and University, P´azm´any P´eter s´et´any 1 / A, H–1117 Budapest, Hungarye-mail:
[email protected] e-mail:
[email protected] Institut f¨ur Astronomie und Astrophysik, Abt. Computational Physics, Universit¨at T¨ubingen,Auf der Morgenstelle 10, D–72076 T¨ubingen, Germanye-mail: [email protected]
Received ; accepted
ABSTRACT
Context.
Based on radial velocity measurements it has been found recently that the two giant planets detected around the starHD 73526 are in 2:1 resonance. However, as our numerical integration shows, the derived orbital data for this system result in chaoticbehavior of the giant planets, which is uncommon among the resonant extrasolar planetary systems.
Aims.
We intend to present regular (non-chaotic) orbital solutions for the giant planets in the system HD 73526 and o ff er formationscenarios based on combining planetary migration and sudden perturbative e ff ects such as planet-planet scattering or rapid dispersalof the protoplanetary disk. A comparison with the already studied resonant system HD 128311, exhibiting similar behavior, is alsodone. Methods.
The new sets of orbital solutions have been derived by the Systemic Console ( ). The stabilityof these solutions has been investigated by the Relative Lyapunov indicator, while the migration and scattering e ff ects are studiedby gravitational N-body simulations applying non-conservative forces as well. Additionally, hydrodynamic simulations of embeddedplanets in protoplanetary disks are performed to follow the capture into resonance. Results.
For the system HD 73526 we demonstrate that the observational radial velocity data are consistent with a coplanar planetarysystem engaged in a stable 2:1 resonance exhibiting apsidal corotation. We have shown that, similarly to the system HD 128311, thepresent dynamical state of HD 73526 could be the result of a mixed evolutionary process melting together planetary migration and aperturbative event.
Key words. planets and satellites: formation, celestial mechanics, hydrodynamics, methods: N-body simulations
1. Introduction
Nearly one third of the multiplanet extrasolar planetary systemscontain pairs of giant planets engulfed in mean motion reso-nances. Resonant extrasolar planetary systems have special im-portance when studying the formation of extrasolar planetarysystems as they require a dissipative mechanism operating on theplanets which is capable of changing the semi-major axis of theirorbits. Hence, resonant planetary systems may support the plan-etary migration scenario, and thus help in answering the ques-tion why Jupiter-sized giant planets do not orbit at their forma-tion place. According to the classical formation theories such asplanetesimal accretion followed by core instability or the grav-itational instability model for the giant planets, it is generallyaccepted that they formed quite far from their host star (beyondthe snowline, which is a ∼ − a = ffi ciently Send o ff print requests to : Zs. S´andor slow (adiabatic), the resonant configuration is preserved dur-ing the migration, and the two planets can travel very close totheir host star. The very e ffi ciency of this mechanism is clearlydemonstrated by the well known resonant system around the starGJ 876. The giant planets in this system are engaged in a 2:1resonance, having very well determined elements due to theirshort orbital periods ( P ≈ P ≈
60 days). It is true in gen-eral cases (Beaug´e et al. 2006) and also in the particular case ofGJ 876 (Lee & Peale 2002; Kley et al. 2005), that as a result ofan adiabatic convergent migration process, beside the mean mo-tion resonance the orbits of the giant planets also exhibit apsidalcorotation, (or in other words apsidal resonance), where the os-culating orbital ellipses of the giant planets rotate with the samemean angular velocity.Numerical integrations based on recently observed orbitaldata of the planetary system HD 128311 (Vogt et al. 2005) showthat the giant planets are in a 2:1 mean motion resonance withoutexhibiting apsidal corotation. In a recent study, S´andor & Kley(2006) o ff ered a mixed evolutionary scenario for this systemcombining an adiabatic migration leading the system into a 2:1resonance and a sudden perturbation. This perturbation could beeither a close encounter between one of the giants and a rela-tively small mass ( ∼ M ⊕ ) planet already existing in the sys-tem (which will refer to in the following as planet-planet scatter-ing), or the fast termination of the planetary migration due to asudden dispersal of the protoplanetary disk. Both these perturba-tions can be strong enough to induce relatively large-amplitude Zs. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 oscillations of the eccentricities and the resonant angles. Whileretaining the mean motion resonance, the planet-planet scatter-ing may even break the apsidal corotation between the orbits ofthe giant planets.In this paper we will support the above described mixed evo-lutionary scenario by the detailed analysis of the newly discov-ered resonant system HD 73526. In a recent paper Tinney et al.(2006) reported that the giant planets of HD 73526 are in the2:1 resonance. It has been found that only the resonant angle Θ = λ − λ − ̟ librates, while both Θ = λ − λ − ̟ and ∆ ̟ = ̟ − ̟ circulate. Here, λ , denote the mean lon-gitude of the inner and outer planet, and ̟ , their periastronlongitude. However, according to our numerical integration, therecently published orbital data by Tinney et al. (2006), as shownin Table 1, result in a weakly chaotic and irregular behavior ofthe system.Since chaotic behavior is uncommon among the known reso-nant extrasolar planetary systems, and may not guarantee the sta-bility of the giant planets for the whole lifetime of the system, wehave searched for regular orbital solutions for the giant planets aswell. By using the Systemic Console ( ),we have found new coplanar orbital solutions for the giant plan-ets around HD 73526 exhibiting regular behavior. The regularnature (and thus the long-term stability) of these orbits have beenchecked carefully by numerical integrations and the RelativeLyapunov indicator chaos detection method (S´andor et al. 2004).Numerical integrations based on our newly derived orbital dataexhibit very similar behavior of the eccentricities of the planetsin HD 73526 to those around HD 128311.This paper aims at a complete analysis of the resonant systemHD 73526 including the stability investigation of the new orbitalsolutions and modeling the formation of the resonant systemHD 73526. The paper is structured as follows: first we presentfour sets of new orbital data. Then we investigate the stability ofthese orbital solutions by using the Relative Lyapunov indicatorchaos detection method, and finally, we o ff er possible evolution-ary scenarios for the system, which are similar to the studiedones in the case of HD 128311.
2. Original orbital data and their stability
It has been mentioned also in the case of HD 128311 byVogt et al. (2005), that a kinematic orbit fit for orbital parame-ters in the case of a strongly interacting resonant pair of giantplanets based only on two unperturbed ellipses, often results ina destruction of the system within a few thousand years. This istrue for HD 73526 as well.Thus Tinney et al. (2006) have also carried out a self consis-tent three-body dynamical orbit fit to the observed velocities (aslisted in Table 1). They found this dynamical fit to be stable over1 Myr integration time. Based on these results of Tinney et al.(2006), we have also performed numerical integration of the sys-tem. We have found that the giant planets around HD 73526 areindeed in a 2:1 mean motion resonance having libration in theresonance variable Θ with an amplitude ∼ ◦ .However, the behavior of the eccentricities is irregular, seeFigure 1. This means that the system is chaotic, which ofcourse may not exclude a practical stability of the giant planets.However, based on the behavior of the observed resonant sys-tems to date, we would prefer orbital solutions exhibiting regularorbits showing a clear dynamical stability for the whole life-timeof the system.By using the Relative Lyapunov indicators (RLI,S´andor et al. 2004), we have mapped out the stability properties e cc en t r i c i t i e s time [years] e e Fig. 1.
Irregular behavior of the eccentricities of the giant plan-ets in the system HD 73526 using as initial conditions the dataprovided by the original dynamical orbit fit (Table 1). The upper(lower) curve shows the behavior of the eccentricity of the inner(outer) planet.
Table 1.
Original dynamical orbit fit of HD 73526 planet mass [M J ] a [AU] e M [deg] ̟ [deg]inner 2.9 0.66 0.19 86 203outer 2.5 1.05 0.14 82 13 χ ν rms [m / s] o ff set velocity [m / s]1.57 7.9 -29.96 of the phase space in the close neighborhood of the originalorbital data set, in a similar way done by Bois et al. (2003) in thecase of HD 160691. We have calculated the stability propertiesof the a − a , e − e , M − M , and ̟ − ̟ parameter planes,where a is the semi-major axis, e is the eccentricity, M is themean anomaly and ̟ is the longitude of periastron of one of thegiant planets, while the indices ‘1’ and ‘2’ refer to the inner andthe outer planet, respectively.In Figure 2 we display the stability structures of the parame-ter planes for the semi-major axis and the eccentricities. Duringthe calculation of a particular parameter plane the other orbitaldata have been kept fixed. This implies that when investigatingthe parameter plane of the semi-major axes a − a , the other or-bital data such as the eccentricities ( e , e ), the mean anomalies( M , M ), and the arguments of the periastrons ( ̟ , ̟ ) havebeen kept fixed to their original values (Table 1). On each pa-rameter plane we have marked the stable regions by white, theweakly chaotic regions by grey and the strongly chaotic regionsby black colors. The actual values of the corresponding orbitaldata are also shown by a symbol ‘ + ’ on each parameter plane. Bystudying carefully the di ff erent stability maps (also the not dis-played but calculated cases of the parameter planes ̟ − ̟ and M − M ), we conclude that the orbital data given by Tinney et al.(2006) are embedded in a weakly chaotic region. We stress thatthis does not imply automatically the instability of their fit, how-ever by using these orbital data the system exhibits irregular (orchaotic) behavior and may be destabilized at later times.We have also studied the stability of the system HD 73526with slightly modified initial orbital elements. We have found,for instance, that even a small error of the initial semi-major axis( δ a = .
01) or eccentricity ( δ e = .
01) of the outer planet leadsto the destruction of the system in a few hundred thousand years.We note that the above errors ( δ a and δ e ) are well inside the s. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 3 Fig. 2.
Stability maps calculated by using the original dynamicalorbit fit listed in Table 1. The structure of the parameter planesindicates clearly that the orbital data given by the original fit(marked by ‘ + ’) are in a weakly chaotic region. Here brighterareas refer to more stable and darker ones to unstable regions.error limits of the original orbit fit of Tinney et al. (2006), being ∆ a = ± .
08 and ∆ e = ± .
09. Thus we believe that the use ofinitial orbital data resulting in regular orbits is more convenientwhen investigating the system HD 73526.On the other hand, the possibility that the system showschaotic behavior as a result of its formation scenario cannot beexcluded. It can also happen that a close encounter with a smallmass body results in a chaotic behavior (characterized by the ir-regular variation of the eccentricities for instance), though thesystem itself remains stable for very long time. However, as ournumerical simulations show, we have not found such a behavior;either the system shows regular behavior, or it is destroyed soonafter the scattering event.
Table 2.
New dynamical orbit fits of the giant planets aroundHD 73526
Fit planet mass [M J ] a [AU] e M [deg] ̟ [deg]1 inner 2.42 0.66 0.26 69.8 206.6outer 2.58 1.045 0.16 163.2 265.62 inner 2.62 0.66 0.209 77.7 208.5outer 2.56 1.047 0.194 131 316.23 inner 2.415 0.659 0.26 70.7 202.9outer 2.55 1.045 0.107 170.7 253.74 inner 2.675 0.66 0.209 80.1 207.1outer 2.53 1.048 0.172 122.4 327.6 Table 3.
Properties of the dynamical orbit fits for HD 73526.
Fit No. χ ν rms [m / s] o ff set velocity [m / s]1 1.81 8.36 -46.672 1.60 8.09 -33.763 1.87 8.4 -49.544 1.58 8.04 -32.33
3. New orbital data and their stability
In order to find orbital solution for the system HD 73526 ex-hibiting regular behavior, we have used the Systemic Console( ), and have found some new orbital fitsfor the giant panets. In Table 2 we list four set of them, display-ing di ff erent behavior. In the case of the first three sets of ourorbital data, beside the mean motion resonance, the giant planetsare also in apsidal corotation but with enlarged amplitudes in theresonant angles Θ or ∆ ̟ , while in the fourth case the apsidalcorotation is no more present. In all four cases the eccentricitiesof the giant planets show relatively large oscillations. The be-havior of the eccentricities is very similar to those found in thecase of HD 128311 indicating clearly, that the present behaviorof the system HD 73526 is not likely to be the result of a smoothadiabatic migration scenario alone.In order to investigate the stability of the newly derived or-bital elements (shown in Table 2), we have computed (by usingthe Relative Lyapunov indicator) a series of stability maps, sim-ilar to those shown in Figure 2. We have found that in all casesthe orbital elements are deeply embedded in the stable regions ofthe parameter planes a − a , e − e , ̟ − ̟ , and M − M . InFigure 3 we show the stability structure of the parameter planes a − a and e − e in the case of Fit 1. The semi-major axesand eccentricities corresponding to Fit 1 are marked by the sym-bol ‘ + ’, and they almost lie in the middle of the stability region(white regions of the parameter planes). We note that at first sightthe structure of the parameter planes a − a corresponding tothe original orbital data (shown in Figure 2) and to the orbitaldata of Fit 1 (shown in Figure 3) are quite similar to each other.However, in Figure 2 the stable regions are very narrow, the lightregions (which may indicate stability) are mainly grey and havea fuzzy structure indicating clearly the weakly chaotic charac-ter of the system. In Figure 3 the light regions are white and theyhave very homogeneous structure showing ordered and thereforestable behavior of the system.We have also checked by numerical integration whether theorbital data shown in Table 2 result in a reliable radial veloc-ity curves indicating the motion of the star. Two examples areshown in Figure 4, corresponding to the Fit 1 and Fit 4, re-spectively. We note that the radial velocity curves of Fit 1 andFit 3 are very similar to each other, while the radial velocitycurves of Fit 2 and Fit 4 are very similar to the radial velocity Zs. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526
Fig. 3.
Stability maps calculated by using the orbital data pro-vided by Fit 1 of Table 2. It can be seen from the parameterplanes a − a ( top ) and e − e ( bottom ) that the orbital dataof Fit 1 (marked by ‘ + ’) are embedded well in the stable whiteregion.curve given by Tinney et al. (2006) corresponding to their orig-inal dynamical orbit fit. In the case of Fit 1 and Fit 3 the massof the inner giant planet is (slightly) smaller than the mass of theouter planet, which is also supported by the planetary migrationscenario. According to full hydrodynamical simulations a cav-ity opens between the giant planets, and the inner planet orbitsin a low density gaseous environment. The outer giant planetis still connected to the disk thus it can accrete more material,which may result in larger mass for it than for the inner planet(Kley et al. 2004).The similarities between Fit 1 and Fit 3 or Fit 2 and Fit 4can also be observed in the behavior of the eccentricities of thegiant planets. The eccentricities in the cases of Fit 1 and 3 showmoderate oscillations, while in the cases of Fit 2 and 4 they os-cillate with very large amplitudes. One example for the behaviorof the eccentricities is shown in the top of Figures 5 correspond-ing to Fit 3. The other two (middle and bottom) panels of Figure -250-200-150-100-50 0 50 100 0 500 1000 1500 2000 2500 r ad i a l v e l o c i t y [ m / s ] JD-2451212 [days]-250-200-150-100-50 0 50 100 0 500 1000 1500 2000 2500 r ad i a l v e l o c i t y [ m / s ] JD-2451212 [days]
Fig. 4.
Radial velocity curves of the star obtained by numericalintegration using as initial conditions the orbital data given byFit 1 (top panel) and Fit 4 (bottom banel) of Table 2. It can beseen that the calculated curves fit well to the measured radialvelocity points (black dots).5 show the behavior of the resonant angles Θ and ∆ ̟ , respec-tively, corresponding again to Fit 3.In the following, when investigating the formation of the sys-tem HD 73526, we intend to model a similar behavior providedby Fit 1 and Fit 3, in which case the outer planet is a little bitmore massive than the inner one, and the oscillations in the ec-centricities are more moderate than in the cases of Fit 2 and Fit4.
4. Formation of the system HD 73526 by smoothinward migration
Since the present behavior of the planetary system aroundHD 73526 is very similar to the behavior of the systemHD 128311, we shall investigate in this section the possibilitywhether this behavior could be the result of a mixed formationscenario combining a smooth migration with a sudden pertur-bative e ff ect as suggested by S´andor & Kley (2006) in the caseof HD 128311. We first model the evolution of the system byapplying a slow inward migration. Having estimated the basiccharacteristics of the planetary migration, we intend to study theformation of the system taking into account (i) a fast dispersalof the protoplanetary disk (resulting in the termination of themigration of the outer giant planet) and (ii) planet-planet scatter-ing. s. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 5 e cc en t r i c i t i e s time [years]e e -80-60-40-20 0 20 40 60 80 0 1000 2000 3000 4000 5000 Θ [ deg ] time [yrs]-150-100-50 0 50 100 150 0 1000 2000 3000 4000 5000 ∆ ϖ [ deg ] time [yrs] Fig. 5.
The behavior of the eccentricities e and e ( top ), the res-onant angles Θ ( middle ) and ∆ ̟ ( bottom ) of the giant planetsby using the initial orbital data given by Fit 3. To study the general feasibility of forming a resonant systemsuch as HD 73526 through a migration process of two embed-ded planets we performed full hydrodynamic evolutions wherethe accretion disk is treated in a flat two-dimensional approxi-mation. The general setup and the applied numerical methods ofthe models are identical to those used in the detailed study ofGJ 876 by Kley et al. (2005), see also Kley (1999). Here we re-strict ourselves to a description of the particulars of our modelfor HD 73526, and concentrate of new features.For the masses of the two planets we use m = . M J and m = . M J (see Fit 3 from Table 2), where M J isthe mass of Jupiter, and the planets are not allowed to accreteany mass during their evolution. The planetary orbits are ini- tially circular at distances of r = r = . M ⊙ . The flat disk extends from of r min = . r max = . Σ ( r ) = Σ r − / . Here, the density at r = Σ r / M ⊙ = . · − which gives a total diskmass in the domain of about 0 . M ⊙ . To make the initial evo-lution of the planets more realistic we superimpose initial den-sity gaps to the unperturbed r − / -profile (Kley 1999, 2000). Thedisk is driven by an α -type viscosity with α = .
01, and the tem-perature is specified using a fixed relative vertical thickness of H / r = .
05, no energy equation is solved. The smoothing lengthfor the gravitational potential is ǫ = . H , which gives approx-imately the correct values of the planetary migration rate whencompared to the full three-dimensional case. For the torque cut-o ff , i.e. the region around the planet which is excluded in the cal-culation of the torques, we use r torq = . R Roche . In the regular(non-resonant) evolution of the planets the exact values of ǫ and r torq should not influence the results too strongly due to the lackof material in the extended gap regions surrounding the planets.When the eccentricities of the planets increase during the evo-lution they periodically move through the disk at their peri- andapo-center and these numerical parameter may become more im-portant. We have checked that our results do not depended on theexact values of ǫ and r torq as long as they are in the above speci-fied range.The computational domain from r min to r max is covered by240 radial and 504 azimuthal gridcells where the radial spac-ing is logarithmic and the azimuthal equidistant. This givesroughly squared gridcells throughout the whole domain. At theouter boundary we use damping boundary conditions where thedisk profiles are relaxed to the initial profile (de Val-Borro et al.2006). At the inner boundary we use a new type of boundarycondition where we specify a outflow condition with an az-imuthally averaged radial velocity which has the magnitude of v r ( r min ) = − ν r min , (1)The typical viscous radial inflow velocity of accretion disks is v visc ( r ) = ν/ (2 r ). We use here a 5 times higher value than v visc ( r min ) to accomodate an increased clearing of the inner diskdue to i ) the disturbances of the embedded planets and ii ) thevicinity of the central star. Nevertheless, this boundary conditionwill prevent material leaving the inner disk too rapidly which istypically observed in hydrodynamic evolutions with embeddedplanet (Crida et al. 2007).In Fig. 6 we show the distribution of the surface density Σ ( r , ϕ ) of the protoplanetary accretion disk after about 1400years. The inner planet is located at x = − . , y = − .
61 andthe outer one at x = − . , y = − .
37. In contrast to previousmodels where the two planets orbit in an inner cavity of the diskwithout any inner disk (Kley et al. 2004), here the inner disk hasnot been cleared and is still present due to the more realistic in-ner boundary condition given in Eq. (1).As a consequence the inner planet still is in contact with theinner disk which exerts gravitational torques on it. These tendto push the planet outward and will produce a damping of theeccentricities. In Fig. 7 we display the evolution of eccentricitiesand semi-major axis of the two planets. Initially, both planets or-bit in a wide joint gap where the outer planet experiences onlythe torques of the outer disk and migrates inward, while the in-ner disk pushes the inner planet slightly outwards. Capture inthe 2:1 mean motion resonance occurs at 400 yrs after whichthe planets migrate jointly towards the star. The eccentricitiesincrease strongly and reach equilibrium values at approximately
Zs. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526
Fig. 6.
Two dimensional density distribution of a protoplane-tary accretion disk with two embedded planets after about 1400years. From their initial positions ( a = . a = . a = .
85 and a = .
36 AU with ec-centricities e = .
32 and e = .
08. Within the inner planetthe inner disk has not been cleared even after such a long evo-lution time due to a modified new inner boundary condition forthe radial inflow velocity.1100 yrs. In contrast additional full hydrodynamical test modelsthat do not have an inner disk show a continued eccentricity in-crease in particular for the outer planet. We attribute this reducedmagnitudes of the eccentricities to the damping action of the in-ner disk. When the eccentricity of the inner planet increases itwill periodically have to enter into the inner disk and the gravi-tational torques will tend to reduce its eccentricity. Equilibriumis then given by the eccentricity driving caused by the inwardmotion of the planets and the damping caused by the disk. Afterresonant capture the resonant angles Θ and ∆ ̟ are both librat-ing with a libration amplitude of ≈ ◦ for Θ and < ◦ for ∆ ̟ . Hence, this hydrodynamical evolution displays an adiabaticmigration process which results in apsidal corotation of the twoplanets, where the apsidal lines are always aligned.In the long term the eccentricity of the outer planet ( e ) re-mains at a level of 0.05 to 0.07 and does not decline any fur-ther. Additional test computations using di ff erent resolutions andboundary conditions confirm this trend. The shown hydrody-namical simulation represents an explorative study to demon-strate that such a damping mechanism of the inner disk may haveoperated. Using this mechanism it might be possible to accomo-date a larger radial migration for the resonant planets than thatfound previously for the case of GJ 876 (Kley et al. 2005) be-cause eccentricties will not be pumped up to very large values.This topic will be investigated in more detail in future work. Hydrodynamical calculations modeling the migration of plan-ets embedded in a gaseous protoplanetary disk require relativelylarge computational time. On the other hand, the damping e ff ect S e m i - M a j o r A x i s E cc en t r i c i t y Time [Yrs]a1 e1a2 e2
Fig. 7.
The evolution of the semi-major axis and eccentricitiesof the two embedded planets for the full hydrodynamical model.The inner planet ( a , e ) is indicated by the darker curves. After400 yrs the outer planet captures the inner on in a 2:1 mean mo-tion resonance and they are coupled during their inward migra-tion. The eccentricities increase rapidly after resonant captureand settle to equilibrium values for larger times.of a protoplanetary disk to the orbital evolution of the embeddedgiant planets can also be modeled conveniently in the frameworkof the gravitational N-body problem by using well chosen non-conservative forces. These additional forces can be parametrizedby the migration rate ˙ a / a and the eccentricity damping rate ˙ e / e ,or by the corresponding e -folding times τ a and τ e of the semi-major axis and eccentricity of the outer planet, respectively (seeLee & Peale 2002; Beaug´e et al. 2006). The relations betweenthe damping rates and e -folding times are ˙ a / a = − /τ a , and sim-ilarly for the eccentricities ˙ e / e = − /τ e . We define the ratio be-tween the e -folding times K = τ a /τ e , (or ˙ e / e = − K | ˙ a / a | ), whichaccording to Lee & Peale (2002), determines the final state ofthe system in the case of a su ffi ciently slow migration.We recall that in the migration scenarios used up to nowmainly the damping of the outer giant planet has been takeninto account. As we have mentioned previously, we have alsosupposed that beside the outer disk there is also an inner disk,and the planets migrates in a cavity between these disks. Thepresence of the inner disk accelerates the inner giant planet,thus forces it to migrate outward and also damps its eccentric-ity. This e ff ect can also be modeled by using a (repelling) non-conservative force also parametrized by ˙ a / a (having positivesign for outward migration), and K or ˙ e / e .In what follows, we present our results obtained in model-ing the formation of the system HD 73526 by a slow migrationprocess. We have studied three cases. In the first case, besidethe damping e ff ect of the outer disk, we have taken into accountthe e ff ect of an inner disk. In the second case we have assumedthat only the outer planet is a ff ected by the outer disk and forcedto migrate inward. Finally, in the third case we have assumedthat after the resonant capture the protoplanetary disk dissapearsgradually inhibiting the further increase of the inner planet’s ec-centricity. In this latter case only the outer planet has been forcedto migrate inward. Our aim is to provide an evolution of the sys-tem which is in accordance with the eccentricity limits given byFits 1-4. s. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 7 s e m i - m a j o r a x e s [ A U ] e cc en t r i c i t i e s time [years] e e a a s e m i - m a j o r a x i s [ A U ] e cc en t r i c i t i e s time [years]a a e e K=10 K=20 K=10K=20 K=15 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 5000 10000 15000 2000000.10.20.3 s e m i - m a j o r a x i s [ A U ] e cc en t r i c i t i e s time [years]a a e e Fig. 8.
Evolution of the semi-major axes and the eccentricitiesof the giant planets during di ff erent migration scenarios: (i) theactions of an outer and an inner disk are modeled, τ a = years, τ a = − × years and K = K =
10 (
Top ); (ii) only thedamping e ff ect of an outer disk is modeled, τ a = years and K =
10, 15, and 20 (
Middle ); and (iii) the resonant capture takesplace right before the disk’s dispersal obeying a linear reductionlaw between 9 × and 1 . × years, τ a = × years,and K =
10 (
Bottom ).In order to compare the results of dissipative N-body cal-culations to those obtained by hydrodynamical simulations, atthe beginning of the migration the giant planets move on cir-cular orbits at distances a = a = M ∗ = . M ⊙ mass star. The masses of the giant planets havebeen fixed to m = . M J and m = . M J , where M J isthe mass of Jupiter. For the inward migration of the outer giantplanet we have applied τ a = years, while for the outwardmigration of the inner giant planet τ a = − × years. The -150-100-50 0 50 100 150 0 5000 10000 15000 20000 r e s onan t ang l e s [ deg ] time [years] Θ ∆ ϖ -150-100-50 0 50 100 150 0 5000 10000 15000 20000 r e s onan t ang l e s [ deg ] time [years] Θ ∆ ϖ -150-100-50 0 50 100 150 0 5000 10000 15000 20000 r e s onan t ang l e s [ deg ] time [years] ∆ ϖΘ Fig. 9.
Evolution of the resonant angles Θ and ∆ ̟ during themigration scenarios corresponding those displayed in Figure 8.ratio between the e-folding times in both cases was K = θ and ∆ ̟ deviate slightly from 0 ◦ and tend finally to 0 ◦ . N-body calculations also confirm the results of the hydrodynamicalsimulations that an inner disk can damp e ffi ciently the eccentric-ity of an inner planet during a permanent migration.Hydrodynamical simulations for HD 73526 also show thatthe lack of the inner disk results in a continuous increase of theeccentricity of the inner giant planet. In what follows, we willstudy which value of ratio K is necessary to stop the increase ofthe inner planet’s eccentricity when no inner disk is considered.(It is known from the hydrodynamical simulations that K liestypically between 1 − Zs. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526
We have applied a damping of the semi-major axis of theouter planet with an e -folding time τ a = year using di ff erent K values ( K =
10, 15, and 20). The time evolution of the semi-major axes and the eccentricities are shown in the middle panelof Figures 8, while the behavior of the resonant angles θ and ∆ ̟ (with K =
10) is displayed also in the middle panel of Figures 9.Studying the behavior of the eccentricities for the di ff erent ratios K one can see that during the migration process the eccentricityof the outer planet is damped su ffi ciently, while the eccentricityof the inner planet is slightly increasing. If the migration is ter-minated when the semi-major axes reach their actual values (at t ≈ . × year), we need K ∼
15 not to exceed the observedeccentricity of the inner planet. (The upper limit for the eccen-tricity of the inner planet is around e ≈ .
3, as it has beenshown by our numerical integrations based on the orbital data ofFits 1-4.) On the other hand, from hydrodynamical simulationswe know that K ∼ −
10, thus although the above values K = K value can be reduced to K = K =
10 by the presence ofan inner disk as shown above by the hydrodynamic simulations.The evolution of the resonant angles Θ and ∆ ̟ (for K = Θ and ∆ ̟ deviate from 0 ◦ , and having reached theirextrema ( ∼ ◦ for Θ and ∼ ◦ for ∆ ̟ ), they tend to 0 ◦ veryslowly. The reason of this behavior might be that the mass ra-tio between the giant planets ( m / m = . m / m ) crit = .
95 given by Lee (2004).If m / m > ( m / m ) crit , the resonant system evolves temporarilythrough an asymmetric apsidal configuration. We note that thedeviation of the angles Θ and ∆ ̟ is the largest in the the abovecase when only the outer giant planet is damped by an outer disk.It is also clear that the migration of the giant planets shouldterminate when they reach their actual positions. If the termina-tion of the migration is a slow process the system reaches a statevery close to a periodic solution of the corresponding three-bodyproblem studied by Psychoyos & Hadjidemetriou (2005), for in-stance.In order to avoid the use of the relatively high values for K ∼ −
20, we have also studied the possibility when the reso-nant capture between the giant planets takes place just before thedispersal of the protoplanetary disk. We recall that this scenario(proposed by Kley et al. 2005) might have occurred in the caseof GJ 876, if there is not acting a strong damping mechanismon the eccentricities. We have supposed that prior to the reso-nant capture the giant planets migrated separately to their orbitshaving a = .
75 AU, and a = . τ a = × yr with K =
10, and captures the inner planet intothe 2 : 1 resonance, see the bottom panel of Figure 8. The reso-nant capture takes place around t ≈ × years, the eccentricityof the inner planet grows rapidly, but before exceeding the limit e = .
3, the protoplanetary disk has already disappeared andthus the migration is terminated. The disk dispersal happens be-tween 9 × − . × years, obeying a linear reduction law.During this simulation the final value of the inner planet’s eccen-tricity remained around the critical value e ≤ .
3, however theratio of the e -folding times was only K =
5. Formation of the system HD 73526 by mixedevolutionary scenarios
In the previous section we have presented three possible migra-tion scenarios that can push the giant planets deep into the 2 : 1resonance not contradicting to the observed upper limits of the e cc en t r i c i t i e s time [years]e e -120-60 0 60 120 0 2000 4000 6000 8000 10000 Θ [ deg ] time [years]-180-120-60 0 60 120 180 0 2000 4000 6000 8000 10000 ∆ ϖ [ deg ] time [years] Fig. 10. Top : Behavior of the eccentricities of the giant planetsduring an inward migration followed a sudden dispersal of theprotoplanetary disk, τ a = × and K =
15. The giant plan-ets orbit initially at a = a = Middle and
Bottom : Behavior of the reso-nant angles Θ and ∆ ̟ corresponding the above scenario. Afterthe fast termination of the migration the resonant angles also be-have very similar to the case displayed in Figure 4.eccentricities ( e ∼ . e ∼ . s. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 9 ios melting together inward migration and sudden perturbativeevents mentioned in (S´andor & Kley 2006). Recent
Spitzer observations of young stars show that the innerpart of the protoplanetary disk may be emptied due to pho-toevaporation induced by the central star (see D’Alessio et al.2005; Calvet et al. 2005). Thus when approaching the inner rimof such a disk, the inward migration can be stopped rapidly(Masset et al. 2006).If the inner rim is located at a given radius of the disk (e.g.it is stationary) the migration cannot be maintained because inthis case the the outer disk, which is the driving agent of themigration, cannot follow the planet. on its inward drift.Based on our experience in modelling the formation of thesystem HD 128311, we have also investigated the above casewhen the termination of the migration (due to the presence ofa stationary inner rim) happens in a very short timescale, here ∆ t d is =
50 years. In this case we apply a faster migration τ a = × years and K =
15. The behavior of the eccentricities andthe resonant angles Θ and ∆ ̟ is shown in Figure 10. It is clearthat this fast termination of the migration can result in a behaviorof the eccentricities and the resonant angles, which is similar tothe observed case shown in Figure 5.The reason of this behaviour is that in the beginning of themigration process, the resonant angles deviate from 0 ◦ , reachtheir extrema, and tend very slowly to 0 ◦ (see also the middlepanel of Figures 9). If the termination of the migration occurswhen the resonant angles still deviate from 0 ◦ , they begin to os-cillate (with an amplitude equal to their values at the end of thedisk’s dispersal) around their equilibrium (which is in this case0 ◦ ). This scenario works also in the cases of slower migrationrates, however in the case of a faster migration the deviation ofthe resonant angles from 0 ◦ is larger, and hence the amplitude ofthe oscillations in e , is noticeably higher.We should note that the above fast migration of protoplan-ets may not be likely in protoplanetary disks. In thick disks fastmigration may occur, but in this case is not easy to find relevantphysical processes leading to the sudden termination of migra-tion. However, the outcome of the sudden stop of migration char-acterized with the above parameters yields results being in goodagreement with the behavior of the system by using as initialconditions the data of Fit 3. The behavior of the eccentricities of the giant planets aroundHD 73526 is very similar to that observed in the case ofthe system HD 128311 and υ Andromedae. In the case of υ Andromedae, Ford et al. (2005) proposed that such a behav-ior is most likely the result of a planet-planet scattering event.Studying the system HD 128311, S´andor & Kley (2006) havedemonstrated that the observed behavior of that system may bethe consequence of a planet-planet scattering as well. In theirrecent paper Tinney et al. (2006) have also suggested that thedynamical behavior of the system around HD 73526, might bethe result of a dynamical scattering event. In what follows, weshall investigate whether the present behavior of HD 73526 canbe modeled by such an e ff ect.We have assumed that in addition to the giant planets, beingalready in a 2 : 1 resonance and migrating towards the hostingstar, a small mass planet ( ∼ M ⊕ ) is also orbiting close to the e cc en t r i c i t i e s time [years]e inner e e -180-135-90-45 0 45 90 135 180 0 5000 10000 15000 20000 r e s onan t ang l e s [ deg ] time [years] ∆ ϖθ Fig. 11. Top : The evolution of the eccentricities of the giant andthe small mass planets during a planetary migration when theinner giant planet captures the small mass inner planet into aprotective 3 : 1 resonance. This resonant configuration inhibitsthe increase of the small mass inner planet’s eccentricity thusthe system remains stable under the whole migration process.
Bottom : The behavior of the resonant angles θ and ∆ ̟ showsthat beside the 3 : 1 resonance, the orbits of the inner giant andthe small mass planet also are in an asymmetric apsidal corota-tion.hosting star in a quasi circular orbit. (This assumption could berealistic as the discovery of a 7 . M ⊕ planet around GJ 876 shows(Rivera et al. 2005)) Initially, the giant planets are far enoughfrom the small mass planet thus they do not influence signifi-cantly its motion. However, as the giant planets migrate inwardapproaching their present positions, they perturb the motion ofthe small planet. As our numerical experiments show at the cor-responding ratio of the semi-major axes, the inner giant planetcaptures the small mass planet into a 3 : 1 resonance. Once thecapture into the 3 : 1 resonance happens, two di ff erent scenarioshave been detected: (i) the eccentricity of the small mass planetgrows initially, then oscillates around the eccentricity of the in-ner giant planet, (ii) the eccentricity of the small mass planetgrows during the whole migration process, which can result inasurvived possible close encounter between the inner small massand the giant planet.In the case (i) the resonant angles of the 3 : 1 resonance θ = λ − λ inn − ̟ − ̟ inn , θ = λ − λ inn − ̟ , and θ = λ − λ inn − ̟ inn librate around ∼ ◦ , ∼ ◦ , and ∼ ◦ , respectively, while ∆ ̟ ∗ = ̟ − ̟ inn also librates around ∼ − ◦ showing that the orbits of the inner giant and of thesmall mass planet are in an asymmetric apsidal corotation. (Herethe index ‘1’ refers to the inner giant planet, while the index ‘ inn ’ e cc en t r i c i t i e s time [years]e inner e e -180-135-90-45 0 45 90 135 180 0 10000 20000 30000 r e s onan t ang l e θ [ deg ] time [years] Fig. 12. Top : The evolution of the eccentricities of the giant andthe small mass planets during a planetary migration when thesmall mass planet is also captured into a 3 : 1 resonance bythe inner giant planet. In this case the resonant configuration isno more protecting, the eccentricity of the small mass planetsincreases during the migration and its orbit can cross the orbitof the inner giant planet.
Bottom : The behavior of the resonantangles θ during the migration process. We note that the otherresonant angles (of the 3 : 1 resonance) are circulating, thus theorbits of the inner giant and the small mass inner planet are notin apsidal corotation.to the small mass inner planet.) This resonant configuration pro-tects the small mass inner planet from the close encounter (mak-ing impossible the planet-planet scattering), and stabilizes its or-bit during the whole migration. The behavior of the eccentricitiesof the planets and the resonant angle θ (for the 3 : 1 resonance), ∆ ̟ ∗ are shown in Figure 11.In the case (ii) the inner giant planet also captures the smallmass planet into a 3 : 1 resonance, however, in this case only θ librates around 0 ◦ (see the bottom panel of Figure 12), theother resonant angles are circulating, thus the orbits of the innergiant and the small mass planets are not in apsidal corotation.During the migration, the eccentricity of the small mass innerplanet grows (shown in the top panel of Figure 12) resulting in amore and more elongated orbits. Thus eventually the orbit of thesmall mass planet will cross the orbit of the inner giant planetand between the planets a close encounter can happen. After theclose encounter the small planet is either ejected from the systemor pushed into a distant orbit from the hosting star.Having performed several numerical simulations, we havefound that the outcome of the resonant capture described in thecases (i) and (ii), depends mainly on the initial positions of the inner giant and small mass planet, and it is not very sensitiveon the mass of the small inner planet (ranging between 0 .
03 and0 . M Jup ). It is very probably that prior to its secular increase, e inn goes through a sudden change in the family of corotationdescribed and analysed recently also for the 3 : 1 resonance byMichtchenko et al. (2006). Thus a detailed analysis would cer-tainly be useful to study the mechanism and the probability ofthe resonant capture into the 3 : 1 resonance leading to eithera protective behavior or a close encounter with the inner giantplanet and a possible ejection of the small mass planet. However,at the present stage of our investigations we restrict ourselvesonly to demonstrate the mechanism that leads to the increase ofthe eccentricity of the small mass inner planet and to a planet-planet scattering event.In our numerical experiments we varied the mass of the smallplanet between 0 .
03 and 0 . M Jup and the migration timescaleas well. We have performed several simulations, and have foundthat once scattering happens, in the majority of the cases thesystem remains captured into the 2 : 1 resonance. However, insome cases we have observed the complete destruction of theresonance, or even escape of the giant planets. We note that ifthe scattering takes place well before the disk dispersal process,the large oscillations in the eccentricities may be smoothed outdue to the damping acting of the disk acting on the eccentric-ity of the outer giant planet. If the scattering happens during orafter the disk dispersal, the oscillations of the eccentricities arepreserved, and the system reaches its presently observed state.We have found that in order to observe the continuous increaseof the small mass inner planet’s eccentricity the migration of thegiant planets should take su ffi ciently long time.Finally, the question may arise whether it is necessary toassume a kind of temporal synchronization between the termi-nation of the migration and the scattering event. We have per-formed several numerical investigations and found in many casesthat the inward migrating giant planets represent a perturbationstrong enough to make the behaviour of the small mass innerplanet unstable. However, the small planet survived for longertimes, and the close encounter happened when the migration wasalready terminated. On the other hand there can also be a natu-ral temporal synchronization. The migration of the giant planetshould terminate when they reach their presently observed posi-tion. If a small mass planet orbits at a distance from the centralstar, which is needed for the resonant capture, its orbit will becaptured into a 3 : 1 resonance, and therefore its eccentricitypumped up making possible a close encounter with the inner gi-ant planet, and finally, its ejection from the system.In Figure 13, which is the magnification of Figure 12, weshow the typical behavior of the eccentricities of the giant plan-ets, and ∆ ̟ during a migration process before and after a planet-planet scattering. The giant planets are started initially on distantorbits having initial semi-major axes a = a = e -folding time of the outer giant’s semi-major axis is τ a = years and K = M inn = . M Jup ,and it orbits at a inn = . t sc = × years.We applied also a linear dispersal process of the protoplanetarydisk taking place between 2 . × and 3 . × years. Thus theoscillations in the eccentricities induced by the scattering eventare not damped out anymore. Since the giant planets are still ina (protective) 2:1 resonance their stability is preserved for thewhole life-time of the system. We note that in contrary to thecase HD 128311, during modeling the formation of HD 73526we have not observed the breaking of the apsidal corotation. It s. S´andor, W. Kley, and P. Klagyivik: Stability and Formation of the Resonant System HD 73526 11 e cc en t r i c i t i e s time [years]e e -180-135-90-45 0 45 90 135 180 26000 30000 34000 r e s onan t ang l e ∆ ϖ [ deg ] time [years] Fig. 13. Top : The behavior of the eccentricities of the giantplanets during a planetary migration process before and after aplanet-planet scattering event taking place around t sc ≈ × years. Bottom : The behavior of ∆ ̟ . The behavior of the eccen-tricities and ∆ ̟ is very similar of those of Fit 3.is highly probable that due to the large planetary masses andrelatively small semi-major axes, the system is more fragile thanHD 128311, and a su ffi ciently strong perturbation which couldlead to the breaking of the apsidal corotation may destroy thewhole system. However, after the scattering event the resonantangle ∆ ̟ shows large oscillations around 0 ◦ .We note that in the above simulations (gravitational three-body problem with dissipative forces) the small mass innerplanet has not been a ff ected by the inner disk. In this case weused a model without an inner disk, which may also be realisticin some cases. However the damping ratio K =
15 is a bit high;smaller K would result in higher eccentricity of the inner giantplanet. As we have already demonstrated in 4.1, the presence ofan inner disk can damp e ff ectively the eccentricity of the innergiant planet. In order to investigate the e ff ect of the inner disk onthe small mass planet we have also performed full hydrodynam-ical simulations using the same initial parameters as describedin 4.1, but in addition to the giant planets, a M inn = M ⊕ innerplanet has also been embedded in the inner disk and started from0 . t =
400 orbital pe-riods (of the inner giant planet) the two massive planets captureeach other in a 2 : 1 resonance slowing down the inward mi- s e m i - m a j o r a x i s [ A U ] time [orbital periods]a inner a a e cc en t r i c i t y time [orbital periods]e inn e e Fig. 14.
Full hydrodynamical evolution of two embedded mas-sive planets and an small additional small mass inner planet.
Top : The behavior of the semimajor axis of the three planets.
Bottom : The behavior of the eccentricities.gration which nevertheless remains faster than that of the innerplanet which is consequently captured at a later time. We shouldalso note that due to the very long computational time, we havenot followed the complete evolutionary way of the small massinner planet. On the other hand it is very probable, that the in-crease of its eccentricity will result in a scattering with the innergiant planet, and finally, ejection from the system. Similarly tothe eccentricity damping mechanism, the final behavior of theinner small mass planet will also be the subject of a more de-tailed analysis investigating inner disks with di ff erent physicalproperties. However, our results to date clearly show that evenin the presence of an inner disk the eccentricity excitation mech-anism of the small mass inner planet (through resonant capture)works well, and enables scattering between the small mass andthe inner giant planet.
6. Conclusion
In this paper we have investigated di ff erent evolutionary scenar-ios, which may have led to the observed behavior of the reso-nant system HD 73526. According to our numerical integration,the original orbit fit provided by Tinney et al. (2006) results inweakly chaotic motion. Thus, we have derived four sets of or-bital elements, which fit very well to the observed radial velocitymeasurements, and also provide, through a protective 2:1 reso-nance, stable motion for the giant planets during the whole life-time of the system. Two of these fits yield very large variationsof the eccentricities of the giant planets, in one case the giantplanets are not even in apsidal resonance. In the case of the other two fits the variations of the eccentricities are more moderate,and beside the mean motion resonance, the giant planets are inapsidal corotation.Since the behavior of the giant planet’s eccentricities in thesystem HD 128311 is very similar to those in HD 73526, wehave investigated whether this system may also be the result ofa mixed formation scenario incorporating a slow migration fol-lowed by a perturbative event. Similarly to the formation of thesystem HD 128311, we have assumed that the two giant planetshave been formed far from their hosting star, and migrated in-ward due to a planet-disk gravitational interaction. During thismigration the planets have been locked in a protective 2:1 reso-nance.Performing full hydrodynamical simulations of two embed-ded giant planets in a protoplanetary disk we have found that thepresence of an inner disk can damp e ffi ciently the eccentricityof the inner giant planet resulting in eccentricity values being ingood agreement with the observations. This mechanism has notbeen investigated yet, and may o ff er a real possibility to avoidthe use of very high values of the ratio K , for example in mod-eling the resonant system around GJ 876. We have additionallyfound that the e ff ect of an inner disk can be modeled very conve-niently by using non-conservative forces. We have investigatedthree convergent migration scenarios that can lead the system tothe deep 2 : 1 resonance, which are all physically realistic.Similarly to HD 128311 we have assumed that the systemHD 73526 su ff ered a sudden perturbation near the end of themigration of the giant planets. This perturbation could either bethe quick termination of the migration, which may possibly beinduced by an inner rim of the disk and an empty region in-side of it, as indicated by some observations of young stars, ora planet-planet scattering event. We have analyzed the results ofan encounter of one of the giant planets with a relatively small ∼ M ⊕ planet, which orbits around the central star at a closerdistance ( a < ff ect of the disk’s sudden dispersal, wehave found that only a very short disk dispersal time ∆ t d isp = ∆ ̟ , ◦ . At momentof the sudden termination of the migration ∆ ̟ ≈ ◦ , and ∆ ̟ will continue to librate with this amplitude around its equlibrium(0 ◦ ). If Θ is also di ff erent from 0 ◦ , it will also librate around 0 ◦ after the scattering event. We note however that in order to ob-tain the present behavior of the system we need relatively fastmigration of the giant planets τ a = × years and K = Acknowledgements.