© ESO, 2012

1. Introduction

Over the past decade the Spitzer Space Telescope (Werner et al. 2004) has provided complementary information about known debris disks in multiple-star systems. Observations of pairs of main sequence stars made by Spitzer revealed stable circumbinary disks around pairs of stars with separations of 0.04–5.31 AU in 14 systems (Trilling et al. 2007, see Table 4 of that paper). In addition, two of these binary systems have an additional companion.

Another example of a debris disk in a multiple star system is the circumbinary disk of HD 98800. This system is a young (<10 Myr) member of the TW Hydra association (Fekel & Bopp 1993). HD 98800 has four stars in two binaries, A and B, travelling around each other in highly inclined and eccentric orbits. According to Furlan et al. (2007), the observations suggests the existence of an optically thick wall at 5.9 AU around the B components and an inner, empty region. However, some optically thin dust orbits the B binary in a ring between 1.5 and 2 AU. In particular, Furlan et al. (2007) have argued that the peculiar structure and apparent lack of gas in the HD 98800 disk suggest that this system is likely already at the debris disk stage.

In later studies of debris disks in triple star systems (Akeson et al. 2007; Verrier & Evans 2008), the quadruple stellar system HD 98800 was studied as a triple system consisting of the B inner binary orbited by a third star. Akeson et al. (2007) investigated the influence of a third star, with an orbital plane initially inclined at 10° with respect to the B inner binary orbital plane, on the HD 98800 disk. The authors found that the inclined stellar orbit could produce a warp in the disk around the B inner binary. Verrier & Evans (2008) considered the third star in three possible current orbital configurations of the HD 98800 system (Tokovinin 1999). The third star is considered to be in an orbital plane initially inclined at  ~143° with respect to the B inner binary orbital plane. The test particles are initially taken on placed in a cloud of orbits. Verrier & Evans (2008) showed gaps in the disk and three distinct stable populations of test particles: a prograde disk, a retrograde disk, and a high-inclination halo.

It is important to study what happens in debris disks to understand the role of stellar companions in planet formation. A stellar companion within 100 AU probably affects the formation of giant planets (e.g., Nelson 2000; Mayer et al. 2005; Thébault et al. 2006; and others).

In this work, we extend the investigations of debris disks performed by Akeson et al. (2007) and Verrier & Evans (2008). We have considered a range of inclinations for the third star and for the particle’s disk. We estimate the stable region in which particles can survive despite the third star’s inclination and eccentricity. Our goal is not to reproduce the current configuration of the debris disk of HD 98800, but to sample the range of possible outcomes for different initial inclinations of the disk and the third star. The astronomical motivation for our study is to investigate whether/where the orbital stability region of a debris disk is possible around the close binary, and to show how the disk structure can be shaped by its dynamical interaction with an eccentric and inclined third star. Investigating the details of these dynamical interactions is beyond the scope of the present work, therefore we will focus the possible impact of the mean motion resonances on the disk because recent studies have mentioned them as a powerful mechanism for instabilities and gaps in the disk, without investigating them in detail (Verrier & Evans 2008; Farago & Laskar 2010). An analytical study of mean motion resonances in this framework is intended to be carried out in the future.

This paper is structured as follows. In Sect. 2, an overview of the empirical expressions of the stability boundary applied here is presented. In Sect. 3, we present the initial conditions for the disk and the triple stellar system assumed here. In Sect. 4, we summarise the relevant results of our numerical simulations and discuss the effects of the inner binary and third star perturbations. In Sect. 5, we present a study on the identification of mean motion resonances in the particle disk. Finally, our conclusions are presented in Sect. 6.

2. Stability boundary

To study the stability boundary of a circumbinary disk in a triple stellar system, we have chosen the HD 98800 system. This system has been modelled as a triple stellar system in which the B binary is orbited by a particle disk and a more distant third star. The A pair is treated as a single star disturber.

The B binary orbit has a semi-major axis of a_B = 0.983 AU and an eccentricity of e_B = 0.7849. The solar masses of the B pair are M₁ = 0.699 and M₂ = 0.582, and the mass ratio μ_B = M₂/(M₁ + M₂) is 0.45 (Boden et al. 2005). The A binary (hereafter, third star) has three possible orbital configurations (see Table 1), and its mass M₃ was assumed to be the sum of the masses of the B pair (M₃ = 1.281 solar masses). The mass ratio μ_A = M₃/(M₁ + M₂ + M₃) is then 0.5 (Tokovinin 1999).

To study the stability of planets in binary star systems, Holman & Wiegert (1999) have derived empirical expressions of stability boundaries for planetary orbits in binary star systems. According to their results, the stability region is delimited by what we called internal (a_I) and outer (a_O) critical semi-major axes. The internal critical semi-major axis is the closest stable orbit to the two stars (or one of the stars) due to the star pair’s circumbinary perturbations. The outer critical semi-major axis is the outermost stable orbit around the star pair (or one of the stars) due to the perturbations of a distant star.

Using the approximated expressions from Holman & Wiegert (1999), we give the values of a_O for each orbital configuration of the third star in Table 2. We found an inner stability boundary of  ~4.04 AU.

Verrier & Evans (2007) have performed numerical simulations on the stability zones for particles orbiting an inner binary in the presence of a more distant third star. They have presented empirical expressions of the stability boundary for coplanar circumbinary stability by modelling the stars as a low-inclination triple system. The empirical expressions for inner and outer borders are written as follows: Table 2 gives the values of a_O for each studied configuration. For the inner stability boundary, we found a value of  ~ 3.87 AU.

Another approximation of the stability boundary for planetary orbits in binary systems can be obtained using the results from Domingos et al. (2006), who numerically studied the orbital stability of a particle in the framework of the elliptic restricted three-body problem. Empirical expressions of the stability boundaries for prograde and retrograde orbits were also derived. The authors have found a_O values that agree well with those from Holman & Wiegert (1999) and Verrier & Evans (2007). The results from Domingos et al. (2006) differ from previous works in that the empirical expression for the outer critical semi-major axis of the particle also depends on the eccentricity e of the particle. In the case of a prograde orbit, the outer semi-major axis is given by (3)where R_H is the radius of the system’s B-A Hill’s sphere. The a_O values for e = 0 for each configuration of the third star are shown in Table 2.

Analysing Eqs. (2) and (3), the a_O value strongly depends on the e_A, a_A and μ_A values when the particle is on circular orbit. This relation may be a result of the fact that the third star’s gravitational perturbation on the particle can change significantly depending on those parameters. This changes the minimum distance between the third star and the inner binary. The third star’s gravitational perturbation on the components of the inner binary can also change, resulting in a perturbation on the particle. Thus, the combination of these two effects must be important to the particle’s dynamical stability. According to Eq. (3), the critical semi-major axis a_O beyond which the particle would not be stable also depends on the e value.

For the inner stability boundary, the expression from Domingos et al. (2006) does not determine the a_I value for the system considered here. In that study, the a_I value represents the semi-major axis of the collision of the particle with the body that it was orbiting. Holman & Wiegert (1999) and Verrier & Evans (2007) have shown that a_I depends on e_B because it corresponds to the semi-major axis related to the centre of mass of the inner binary, from which the orbit of the particle can become unstable. However, the eccentricity of the particle was not considered in those studies.

Considering that particles can be in high-e orbits with a pericentre distance of R_p = a(1 − e) within the unstable region, particles might be removed by perturbations from the inner binary. To estimate the value of e as a function of the minimum semi-major axis a, in which the particle could move near the inner binary, we have assumed that R_p is given by Eq. (1). The eccentricity of the particle as a function of a is then given by (4)we notice that the difference between the a_I values from Holman & Wiegert (1999) and those from Verrier & Evans (2007) were not significant. The advantage of our approach is that the a_I value depends on not only μ_B, a_B and e_B, but also on a whole set of particle eccentricities.

Assuming a limit of eccentricity, e_lim, for an orbit with a semi-major axis limit, a_lim, to be within the stable region, it was possible to establish the eccentricity from expressions (3) and (4). The (a_lim, e_lim) maximum pair was estimated by setting R_p equal to the outer stability semi-major axis, a_O. For the I, II and III configurations, these values are approximately (8.830 AU, 0.572), (6.670 AU, 0.433) and (5.963 AU, 0.366), respectively. For the a = a_lim case with e exceeding e_lim, the pericentre and/or apocentre distances can be located outside the stable region. When e = e_lim, the particle orbits have pericentre and apocentre distances on the inner and outer stability boundaries, respectively.

According to the goal of this work, we focused on particles within the stable region. In our analyses, we used the definition of the (a_lim, e_lim) pair as a reference value. These physical quantities represent the dynamical effects of the inner binary and third star companion on a circumbinary disk and consequently represent the particle stability. The limit radius is between a_I and a_O values. The combination of e_lim and a_lim represents a critical value of high eccentricity and semi-major axis values of orbits relative to the particle instabilities. A specific (a, e) pair can result in an orbit of collision with the inner binary or escape out of the system.

2.1. Critical inclination

It is particularly interesting to determine the conditions under which an inclined protoplanetary disk remains stable. The mean inclination i of a disk might also have implications for the stability region of a system in which large bodies disturb the disk, producing dust through collisional events. Depending on the system characteristics, something will disturb the disk’s increasing i value (Artymowicz 1997; Thébault & Augereau 2007; Quillen et al. 2007). An important question regarding this problem is the existence of a critical value for the inclination ( ~ 39.23°) between the orbital planes of the disturbed and disturbing bodies (restricted three-body problem). If the disturbing body is in a highly inclined orbit relative to the orbital plane of the disturbed body, the eccentricity increases and the near-circular orbit becomes highly elliptic. This effect is called the Kozai-Lidov mechanism (Kozai 1962) and (Lidov 1962).

The Kozai-Lidov mechanism causes strong periodic variations of the particle eccentricity and inclination while conserving their semi-major axes. The concurrent increase in the inclination and decrease in the eccentricity usually leads to instability in the system. Assuming that the particle’s initial eccentricity is close to zero, we can obtain the maximum eccentricity of a particle given an initial inclination (Innanen et al. 1997). Carruba et al. (2002) have shown that particles on higher inclination librating in the Kozai-Lidov resonance may not reach the maximum eccentricity and can be stable for millions of years.

In the system studied here, the Kozai-Lidov mechanism takes effect on the circumbinary disk when the relative inclination of the particle and third star orbital planes is greater than 40°. Verrier & Evans (2009) showed that the inner binary may cause a nodal libration instead of Kozai-Lidov cycles, which stabilises the test particles against any Kozai-Lidov instability driven by the outer star. These authors also reported that particles closer to the inner binary have nodal libration periods shorter than Kozai-Lidov cycles. Therefore, the nodal libration dominates, resulting in the stability of the particle. Otherwise, when particles more distant (closer to the outer border of the stability region) have Kozai-Lidov periods shorter than the nodal libration, they are destabilised and ejected from the system or collide with the inner binary.

3. Numerical simulations

Using the three possible fits for the third star given in Table 1, we numerically investigated the stability of the particles around the B pair. The B binary orbit was assumed to be in the plane. The orbit of the third star was placed around the centre of mass of the B pair. The initial inclinations for the third star (I_A) are given with respect to the initial B binary plane and are taken to be from 0° to 50° with ΔI_A = 10°.

The disk was formed by 30 000 massless particles (planetesimals) distributed on the same plane around the B binary with semi-major axis (a) from 2 to 30 AU. The initial inclination (i) of the disk was assumed to vary from 0° to 90° with Δi = 10° with respect to the initial B binary plane. The initial eccentricities (e) of the particles were randomly distributed within 10^-3. Longitudes of pericentre and mean longitudes were distributed randomly between 0° and 360°. Our reference plane is the B orbital plane. In particular, our main interest is the final outcome of the evolution of the disk with respect to the initial B orbital plane at time t = 0.

The numerical simulations were stopped after 1 Myr and interrupted if any of the following situations occurred: (i) the planetesimal travelled closer than 0.01 AU from B binary components or (ii) the planetesimal travelled farther than 200 AU from B binary (ejection of the system). The particles surviving the full integration time were considered stable. We used the Swift package (Wisdom & Holman 1991; Levison & Duncan 1994) and the Hierarchical Jacobi Symplectic algorithm (Beust 2003).

4. Results and discussion

Our numerical results are summarised in Figs. 1 to 4, which refer to the final orbital evolution of the particles on the a − e and a − i planes. Similar results have been obtained for three orbital configurations. As an example, Figs. 2 to 4 show the final results for the orbital configuration I.

4.1. The coplanar case

Figure 1 shows the final results of the numerical simulations for orbital configurations I, II and III for the coplanar regime (I_A = i = 0). In configuration I, the outer critical semi-major axis is much larger than that in configurations II and III. Figure 1 shows the final eccentricity as a function of the final semi-major axis of the particles. The stable region is delimited by the two borders, a_I and a_O. The left border (red curve) and the right border (green curve) are given by Eqs. (4) and (3). Table 2 gives the values of a_O obtained for the simulations for each orbital configuration. The inner critical semi-major axis at and beyond which particles survive the integration is  ~3.3 AU.

In Fig. 1, it can be noted that particle orbits were not restricted to within a (a, e) triangular region. Orbits within this region are considered to be stable, and orbits located outside this region are considered to be unstable. It is important to note that the simulations were halted after 1 Myr, which corresponds to 10% of the age of a system such as HD 98800. If the simulations were continued, some of the remaining particles whose orbits were outside of the stability region would have also been ejected or collided.

The surviving particles within the stable region (hereafter the “protected region”) are expected not to be ejected or even collide with the central stars. All particles were initiated with near circular orbits; at some point, some particles achieved eccentric orbits. For orbital configuration I, some particles are excited to higher eccentricities around a =  4.0, 4.2, 10.0, 10.5 and 11.5 AU. For orbital configurations II and III, particles in high-e orbits appear around 5.0 AU and 5.5 AU, respectively, and appear close to the outer border. In general, when particles are not in high-e orbits, they remain in orbits with e ≤ 0.1 within the “protected region”.

Our results have indicated high sensitivity to the third star’s eccentricity. For configurations with a high eccentricity of the third star (0.5 and 0.6), particles in the outer region of the disk (a ≥ 8.0 AU) are ejected from the system on a short timescale. For large eccentricities of the third star, the value of e_lim decreases, and the timescale for particles reaching this limit must be shorter than that for systems with low eccentricities of the third star. The three estimates of the value of e_lim are in reasonable agreement with the results of our numerical simulations.

4.2. Stellar perturbations on inclined orbits

In this section, we present results for inclined orbits of the third star and the particle disk. First, we have considered the case of an inclination between the B pair orbit and the particle’s disk, followed by an inclination between the disk and the third star orbit.

4.2.1. Inclined particle disk

Figure 2 shows that for the initial inclinations of 10° and 20° the final results in a − e planes are identical to those shown in Fig. 1. The final distribution of the eccentricity of the particles within the “protected region” also does not significantly differ from those in Fig. 1. The particle concentrations in high-eccentricity orbits have not changed location with respect to the coplanar case. Figure 2 shows the final inclination of the particles. For the disk initially inclined by 20° the distribution of maximum final inclination of the particles was  ~ 60°. Despite a higher final inclination, the final distribution of the number of particles is similar to that for the coplanar case. In particular, for configuration I, we note that the particle concentration for the higher inclination occurs at a = 4.5 AU.

Figure 3 shows the final orbital eccentricities and inclinations of the particles for a disk initially inclined by 30°, 50° and 90°. Evidently, a significant number of particles have their eccentricities and inclinations excited to very high values for the disk initially inclined by 30°. We have seen that many particles were quickly ejected in the vicinity of the outer boundary. In the inclination plot, we can note that between  ~4.0 and 7.5 AU particles are excited in inclination up to  ~160°.

For particle disks with initial inclinations 50° and 90°, results showed that some particles remained stable for high values of initial inclination. Note that these particles have significantly increased their eccentricities, but the eccentricities did not reach sufficiently high values to be removed from the system during the integration time. On the other hand, it can be noticed that many particles have not significantly increased eccentricities.

We can also note that for i = 90° particles are stable up to a ≲ 9.0 AU, with gaps appearing around  ~2.5 AU and 3.1 AU. In general, most of these particles have final eccentricities of e ≤ 0.1 when the disk orbit is inclined. Depending on the initial value of i, the a_I value at and beyond which particles survive the integration is within 2.4 to 3.3 AU. The lower value for a_I was obtained by considering high relative inclinations (i > 50°).

4.2.2. Relative inclination between the particle disk and third star

Another perspective is to initially consider the particle disk and third star on non-coplanar orbits. Mazeh & Shaham (1979) have shown that when the inclination of the third star is nonzero, the eccentricity of the inner binary oscillates (with corresponding oscillations in its inclination), while I_A and e_A remain approximately constant. As expected, our numerical simulations have shown that for higher values of I_A ( ≳ 40°), the effect of the inner binary perturbations on the inner edge of the particle disk becomes important in determining the inner stability boundary.

Figure 4 presents a representative sample of the final results when the particle disk and third star were placed in non-coplanar orbits. The initial I_A value was 30°, and the initial inclinations of the particle disk were 20°, 50° and 90° (indicated in the upper right corner of the plot). Varying i and I_A, the closest distance to the binary at which particles survived was 2.4 AU. As shown in the figure, some results revealed gaps (around 4.5 AU and 5.5 AU) and a diffusion of the particles, with subsequent particle ejection or collision with the inner binary.

In general, we found three distinct and clearly located structures: (i) a chaotic region closer to the inner binary, (ii) a stable region where particles on highly inclined orbits can survive, and (iii) an unstable region that depends on the distance of the third star and the relative inclination of the particles. Depending on initial conditions from the third star and the particle disk, eccentricity peaks and gaps can appear within the stable region.

5. Mean motion resonances

Our investigation was motivated by the eccentricity peaks shown by some particles with low inclination (i ≤ 30°), which might be due to resonant effects. Thus, instabilities and gaps appear around these resonances as a result of high eccentricities.

We are studying a system with four bodies, none of which can be considered a central body. Thus, there are four natural frequencies in the system given by the mean motions of each body. Given that the stars have highly eccentric orbits and that this section is restricted to the planar case, the resonant angle involves not only mean motions but also the pericentre longitude of the particle. The expression of the resonant angle associated with the resonance is then given by (5)where j’s are integers with j₁ + j₂ + j₃ + j₄ + j₅ = 0 according to the first d’Alembert rule. The λ₁, λ₂, λ₃ and λ₄ are mean motion longitudes that correspond to the particle, the inner binary (Ba and Bb stars), and the third star, respectively, and ϖ₁ is the pericentre longitude of the particle. A particle is involved in a mean motion resonance when the critical argument rate corresponds to , where and denote the mean motions of the inner binary, and are the mean motions of the particle and the third star, respectively, and is the pericentre rate of the particle.

It is well known (see for instance Murray & Dermott 1999) that the orbital periods of each star of the B pair about its centre of mass are the same. The mean longitudes in their orbits differ by 180°. Therefore, there is only one frequency (the mean motion) associated with the B pair, rather than two.

Taking , we rewrite in the following form: (6)in Fig. 1, we can see peaks of eccentricity in regions close to the inner binary (3 AU  ≤ a ≤ 5 AU) and in regions farther away (9 AU  ≤ a ≤ 12 AU). Let us call the first region α and the second region β.

In region α, the gravitational effects of the third star do not destabilise the particle. The third star orbital frequency () is much lower than the inner binary orbital frequency (), and thus contributes very little to the combination given by Eq. (6). In practice, we can consider j₄ = 0 such that (6) becomes . This equation involves two-body resonances because we only used two mean motion frequencies, but physically we have three bodies involved in the resonance with λ_Bb = λ_Ba + 180°.

In contrast, the third star plays an important role in region β. These interactions again involve two-body resonances, with the particle seeing the inner binary essentially as a single body. Now we can consider j_B = 0 such that (6) becomes .

Table 3 presents a list of the two-body mean motion resonances obtained here. In the first and second columns, we show the resonant argument and resonant semi-major axis of the particle, respectively. Figure 5 reproduces Fig. 1 with the locations of these two-body mean motion resonances. In Fig. 6, we show examples of temporal evolutions of a, e, and φ for two particles (P1 and P2) identified in the libration regime of the resonant angle.

The particle P1 is an example of a particle quite close to the inner binary. We note that the combination of the resonant angle φ = 8λ₁ − λ_B − 7ϖ₁ (second plot) shows that the particle is in the 8: − 1 two-body mean motion resonance with the inner binary. The time variations of the resonant angle show that it initially librates for approximately 3240 years, at which point the behaviour of φ, a and e clearly indicate the escape of the particle. An important point to note is the manner in which the resonance appears on the plot. There is no obvious pattern to the libration of the resonant angle, but this fact alone does not mean that the orbit is chaotic. The semi-major axis and eccentricity oscillations exhibit good correlation. Note that the variations in a and e are not different. After approximately 1900 years, the evolution of the eccentricity undergoes stronger irregular variations from 0.08 to 0.23 and the value of a remains close to the resonant value with increasingly erratic oscillations up to approximately 3240 years, when the particle escapes the system. Thus, we can say that the dynamical behaviour of this particle is governed by the 8: − 1 two-body mean motion resonance. To verify the importance of the frequency in the resonant angle, we plotted the resonant angle using only and (setting j₅ = 0). We verified that the plots are almost identical. Therefore, the contribution of  to the resonant angle is not significant.

The P2 particle is an example of a particle near the outer border of the stable region. In this resonance, the particle sees the inner binary essentially as a point mass. We show our results for 6000 years to provide a better graphic representation of the resonant angle, but the particle survives the full integration time. The resonant angle is φ =  −λ₁ + 11λ₄ − 10ϖ₁, and the resonant semi-major axis is 9.9331 AU. According to the results, there is no obvious pattern to the libration of the resonant angle. The particles remain oscillating for all times among several values of the resonant angle, and there is no tendency for a circulating regime. The semi-major axis and eccentricity oscillations exhibit a good correlation with each other. To verify the importance of  on the resonant angle of P2, we plotted the resonant angle using only and (setting j₅ = 0). We verified that the plots are quite different. Therefore, the contribution of has a significant effect on the resonant angle and cannot be neglected.

In our study of mean motion resonances, we concentrated on the planar case. The resonant angle was defined using only the mean and pericentre longitudes. Certainly, when the disk and/or the third star is inclined, a study considering the node and pericentre longitudes would be required to gain a better understanding of the resonant effects.

6. Final remarks

In this work, we obtained the borders of the stable region for a particle disk around a close binary system disturbed by a gravitational field of a distant (minimum pericentre distance  ≲ 43 AU) companion star.

Our numerical results showed that the orbital evolution of the particle disk has been affected by the nodal libration due to the inner binary, the Kozai-Lidov effect and several mean motion resonances but is not dominated by any one of these factors in particular. In reality, the situation is much more complex. There are several mechanisms, overlaps of mean motion and/or secular resonances that increase the particle eccentricities, inducing instability. We have shown the following:

• 1.

  There is an eccentricity limit for a particle to remain in thestability region. This limit depends on the a value of the particle. We recall that the value of the particle eccentricity can vary significantly due to the resonant effects (mean motion and secular resonances). For particles at high inclination, the orbits can be destabilised by the Kozai-Lidov mechanism. Depending on the initial conditions, eccentricity and inclination peaks appeared in regions where the effects of the nodal libration influence the stability of the particle orbit. We also showed that this eccentricity limit depends on the system parameters and should be considered when estimating the stability of orbits in triple systems.

• 2.

  Among the mechanisms that can increase particle eccentricity, we have identified two-body mean motion resonances as a powerful mechanism to increase the particle eccentricity on short time-scales (10 kyr or less) and produce a longer-lived high-eccentricity population. These resonances are of high order and are considered a weak class of resonances in the outer asteroid belt; thus, they have not previously been considered in stellar systems. A study of the dynamic behaviour of particles in mean motion resonances is in progress and will be the subject of a future paper.

Although we have used the specific parameters for HD 98800, we hope that our results can contribute for understanding how stellar companions can affect the evolution of disks around close binaries for any system.

Acknowledgments

The authors thank Hervé Beust for providing his HJS algoritm and James Freddy L. Machuca for helpful discussions and comments on resonant dynamics. This work had financial support from UNESP, CNPq and FAPESP (grant number 2008/08679-4). These supports are gratefully acknowledged.

References

All Tables

All Figures