A&A 439, 585-594 (2005)
DOI: 10.1051/0004-6361:20042441

Dynamics of the young multiple system GG Tauri

I. Orbital fits and inner edge of the circumbinary disk of GG Tau A

H. Beust 1 - A. Dutrey 1,2


1 - Laboratoire d'Astrophysique de Grenoble, Université J. Fourier, UMR 5571 du CNRS, BP 53, 38041 Grenoble Cedex 9, France
2 - Laboratoire d'Astrodynamique, d'Astrophysique et d'Aéronomie de Bordeaux, 2 rue de l'Observatoire, BP 89, 33270 Floirac, France

Received 26 November 2004 / Accepted 1 April 2005

Abstract
The GG Tauri system is a prototype example of a young, pre-main sequence multiple stellar system. Both millimeter and optical observations reveal that the central binary (GG Tau A) harbors a ring-shape circumbinary disk. In this paper, we analyse a compilation of astrometric data of the GG Tau A binary over 12 years that allow us, in combination with information coming from the disk observations (millimeter data), to derive a fairly accurate fit of the orbit (a=32.4 AU, e=0.34). We also perform a dynamical study of the circumbinary disk, and show that in order to be in agreement with the orbit deduced from the fit, the inner gap of the disk should be approximatively twice as small as observed. We show that if we allow the error bars on the astrometric data to be larger, another orbit may be found (a=62 AU, e=0.35), compatible with the location of the inner disk gap, although only marginally compatible with the astrometric data. Possible solutions to solve this discrepancy are discussed.

Key words: stars: circumstellar matter - stars: individual: GG Tau - methods: numerical - celestial mechanics - planetary systems: protoplanetary disks

1 Introduction

GG Tauri is one of the best-known multiple T Tauri systems found in the Taurus-Auriga clouds. It is a quadruple system consisting of two binaries. GG Tau A (the first binary) is brighter and closer (0.25 $^{\prime\prime}$). The second pair (GG Tau B) is wider (1.48 $^{\prime\prime}$) and located 10.1 $^{\prime\prime}$ to the south (Leinert et al. 1993). A circumbinary disk orbiting GG Tau A has been spatially resolved in both millimeter and near infrared wavelength domains (Dutrey et al. 1994; Roddier et al. 1996, hereafter R96; Guilloteau et al. 1999, hereafter GDS99). The disk has been extensively studied as a case study of a young binary system. Interferometric observations at 1.3mm of the optically thin thermal dust emission reveal that almost 70% of the material is confined in a sharp edged ring-like structure extending between 180 and 260 AU around GG Tau A, while the 13CO(2-1) line map shows that the rest of the material extends up to 800 AU or more (GDS99). Near-Infrared images performed at the CFHT by R96 have also revealed the circumbinary ring with observed properties (inner radius, position angle and inclination) in excellent agreement with those found by mm interferometry. Table 1 summarizes the observed disk properties.

There is evidence for Keplerian motion of the circumbinary disk around GG Tau A. Kawabe et al. (1993) first noted a velocity gradient which may be indicative of rotation. A detailed analysis of the velocity maps of the disk in 13CO(1-0) and 13CO(2-1) emissions led Dutrey et al. (1994) and GDS99 to conclude that the motion of the disk was essentially Keplerian, and to derive a strong constraint on the central total mass ( $1.28\pm0.07~M_\odot$; GDS99). Meanwhile, White et al. (1999), on the basis of stellar evolution models and photometric data from the HST, gave independent estimates for the masses of the four components of the systems. His determinations for GG Tau A are compatible with the dynamical mass of GDS99. GG Tau B appears to be much less massive than GG Tau A, GG Tau Bb being even a substellar object (0.044 $M_\odot$). Table 2 summarizes the mass estimates of the binary components.

The relative motion of the components of GG Tau A was observed over several years (see Table 3), providing additional kinematic constraints which allowed several authors to derive orbital solutions (R96; McCabe et al. 2002, hereafter MC02; Tamazian et al. 2002).

The information provided by the relative motions of GG Tau A over several years, combined with the knowledge of the geometry and kinematics of the circumbinary disk and the stellar parameters, makes GG Tau A an ideal system for applied studies of the dynamical interaction between a binary and its surrounding circumbinary disk. This problem has been theoretically investigated in many studies, often using smoothed particle hydrodynamics (SPH) (Artymowicz & Lubow 1996,1994; Bate 2000). These studies yielded estimates for the ratio of the radius of the inner gap of the disk to the semi-major axis of the binary. These results were recently confirmed by symplectic integrations using the HJS software (Beust 2003), which are to be compared to SPH integration by Artymowicz & Lubow (1994) with negligible viscosity. The HJS (Hierarchical Jacobi Symplectic; Beust 2003) is a variant of the popular symplectic integration method WHM (Wisdom-Holman Mapping; Levison & Duncan 1994; Wisdom & Holman 1991), but designed for the dynamics of hierarchical stellar N-body systems, while the original method accounts for planetary system dynamics.

Table 1: Observed properties of the GG Tau A binary system.

Therefore, the aim of this paper is to investigate the dynamics of the GG Tau A system (inner binary + circumbinary ring). In Sect. 2, we review all the astrometric data of GG Tau A and derive a new orbital solution that appears to be compatible with the latest data. In Sect. 3, we use the HJS integrator with our orbital solution to investigate the dynamics of circumbinary material orbiting GG Tau A. We show that there is a discrepancy between the orbital fit and the observed location of the inner edge of the disk at 180 AU, unless we allow the error bars on the astrometric measurements to be larger than given in the various references. We discuss this possibility in Sect. 4 and suggest other ways to solve the discrepancy. We present our conclusion in Sect. 5.

We only discuss here the inner profile of the disk surrounding GG Tau A, a probable result of interactions with the inner binary. GDS99 have shown that 70% of the disk mass is located in a narrow ring of width $\Delta R \simeq 80$ AU which also has very sharp edges ($\sim $$10 \pm 5$ AU). The sculpting of the outer edge of the disk is another problem that may be related to dynamical interactions of the GG Tau A system with the second binary GG Tau B. This issue will be investigated in a forthcoming paper (Beust & Dutrey 2005).

Table 2: Masses of GG Tau Aa and Ab used in this paper, plus masses of the GG Tau Ba and Bb components.

2 Orbit solution for GG Tauri A

Table 3: List of astrometric observational data (position angle and separation) for the central binary of GG Tau we use to fit its orbital motion.

The orbital solution of the central binary (GG Tau A) is theoretically well constrained: first, assuming that the circumbinary disk lies in the same plane as the orbit, its inclination angle with respect to the plane of the sky is known ( $37\pm1\hbox{$^\circ$ }$; GDS99); second, its dynamical total mass is constrained by the Keplerian motion of the disk ( $1.28\pm0.07~M_\odot$; GDS99); finally, the astrometric monitoring of the binary performed for many years (see Table 3) provides additional kinematic constraints.


  \begin{figure} \par\includegraphics[angle=-90,width=8.4cm,clip]{2441f1a.ps}\hspace*{2mm} \includegraphics[angle=-90,width=8.4cm,clip]{2441f1b.ps} \end{figure} Figure 1: Least square fits of the projected separation ( left) and of the position angle ( right) of the GG Tau A binary. The data are taken from Table 3. Analysis similar to that of R96.
Open with DEXTER

We know from the disk images that the inclination angle i of the disk with respect to the plane of the sky is $i=37\pm 1\hbox{$^\circ$ }$. The observations reveal that the disk is seen from the south side (GDS99). Following the usual dynamical convention, we will call $\lambda $ the oriented inclination of the disk (or the binary assuming the system is coplanar), which in our case is $\lambda=180\hbox{$^\circ$ }{-}37\hbox{$^\circ$ }\pm 1$; this means that the motion of the disk is seen to be retrograde (see Appendix, Eq. (A.8)). We note $\phi $ the position angle of the small axis of the elliptic shape of the disk on the sky with respect to north (this is the projection of the rotation axis of the whole system onto the plane of the sky when the binary and the disk are coplanar). The data reveal that $\phi=7\pm 2\hbox{$^\circ$ }$ (GDS99).

What is measured on the central binary is its projected separation p, its position angle $\psi $ (with respect to the north), and the temporal derivatives of these quantities. The separation is roughly $0.25\hbox{$^{\prime\prime}$ }$, which leads to p=35 AU assuming a distance of d=140 pc for the Taurus star-forming region (Elias 1978). The position angle $\psi $ has been observed with good accuracy to continuously decrease from $+9\hbox{$^\circ$ }$ to $-14\hbox{$^\circ$ }$ between 1991 and 2002. A compilation of astrometric data gathered over the 13 past years is given in Table 3. The fact that $\psi $ is seen to decrease shows that the disk is viewed from south side (see Appendix), in agreement with the optical and mm data (GDS99).

The binary has been observed around $\psi=0$ in the past years. Once $\lambda $, $\phi $ and M (the binary mass: M=M1 + M2) are known, the knowledge of the other astrometric quantities ( $p, {\rm d}p/{\rm d}t, {\rm d}\psi/{\rm d}t$) allows us to derive the orbital elements a, e and m of the binary, where a is the semi-major axis, e is the eccentricity, and m is the current mean anomaly. As usual in celestial mechanics, the mean anomaly is a quantity proportional to the time that gives the position of the binary on its orbit. We have m=0 at periastron and $m=2\pi$ one orbital period later.

The problem with this technique is that the uncertainties on the fits of the rates of change ${\rm d}\psi /{\rm d}t$ and ${\rm d}p/{\rm d}t$ can lead to very large error bars on the orbital elements. From their data, R96 derive ${\rm d}\psi/{\rm d}t=-2.0~\hbox{$^\circ$ }/\mbox{yr}$and ${\rm d}p/{\rm d}t = -0.4~\mbox{mas}~\mbox{yr}^{-1}\times d$. With $M=1.28~M_\odot$, this yields

 \begin{displaymath}a=111.4~\mbox{AU},\quad e=0.692\;,\quad m=-9.6\hbox{$^\circ$ }. \end{displaymath} (1)

This orbit appears indeed very eccentric, the binary currently being close to periastron. But more recently, adding additional data, MC02 derive e=0.32 and a=35 AU with large error bars, the binary being now closer to apoastron. Independently, using a similar set of data, Tamazian et al. (2002) get e=0.317 and a=36 AU.

Such a discrepancy led us to make a new fit using all the data gathered by various authors over 12 years. These data are listed in Table 3, and the corresponding error bar weighted least square fits are displayed in Fig. 1. We derive

 
$\displaystyle \frac{{\rm d}\psi}{{\rm d}t} = -1.44\pm0.15~\hbox{$^\circ$ }/\mbox{yr}$      
$\displaystyle \frac{{\rm d}p}{{\rm d}t} = 0.28\pm0.76~\mbox{mas}~\mbox{yr}^{-1} \times~d .$     (2)

While the fit of ${\rm d}\psi /{\rm d}t$ can be considered as relevant, obviously  ${\rm d}p/{\rm d}t$ cannot be constrained, as can be seen from Fig. 1. The error bar nevertheless provides upper limits to the variation rate (decrease or increase) of p. With these values and assuming d=140 pc, we derive

 \begin{displaymath}a= 32.4~\mbox{AU} ,\quad e=0.34\;,\quad m = 179.9\hbox{$^\circ$ }, \end{displaymath} (3)

i.e., something much more comparable to the fits of MC02 and Tamazian et al. (2002) than to the one by R96. According to our fit, the binary is currently just before apoastron.
  \begin{figure} \par\mbox{ \includegraphics[width=8cm]{2441f2a}\hfil \includegraphics[angle=-90,origin=br,width=8.8cm]{2441f2b} }\end{figure} Figure 2: Variations of the fitted orbital elements (a,e) of the central binary as a function of the observational parameters. The left plot is a greyscale map of the eccentricity e for d=140 pc, as a function of ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$. The point quoted "R96'' is the R96 value, the "McC02'' point corresponds to MC02 (with its error box, 1$\sigma $), and the "B03'' point corresponds to this work. The right plot shows the zone in (a,e) space corresponding to the error box (1$\sigma $) around the "B03'' point, for three different values of d. The grey shaded area is the total error box in (a,e) space if we let all parameters vary within their error boxes. The two points superimposed (B03, McC02) have the same meaning as in the left plot, and hold for d=140 pc.
Open with DEXTER

However, the orbital elements we derive are very sensitive to the uncertainties on the measurements of the parameters involved in the fit. The orbital elements are functions of p, ${\rm d}\psi /{\rm d}t$, ${\rm d}p/{\rm d}t$, $\lambda $, $\phi $ and M (see Appendix). The major uncertainty on them comes from the uncertainty on the fits of the rates of change ${\rm d}\psi /{\rm d}t$ and ${\rm d}p/{\rm d}t$. But some of the parameters (p, ${\rm d}p/{\rm d}t$, but also M) implicitly also depend on the assumed distance for the Taurus star forming region, which is only known within $\pm $10 pc (Kenyon et al. 1994), at best. This error bar also appears as an additional but significant source of uncertainty on the orbital elements. The sensitivity of the fit to the error bars on the parameters is so high that deriving error bars for the orbital elements from the error bars of the observational parameters is meaningless, as their variation is far from being linear within the error box of the parameters. We thus adopted a non-linear treatment, computing all the possible set of values for the orbital parameters when we let the parameters vary within their error boxes.

This is illustrated in Fig. 2. In the left plot of Fig. 2, we plot a 2D map of the eccentricity of the orbit as a function of ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$, all other parameters being frozen to their mean values (in particular, d=140 pc). In the right plot, we show the error boxes in (a,e) space we derive if we let the parameters vary within their own error boxes. More specifically, in the three outlined boxes, all parameters but ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$ are frozen (each box corresponds to one value of d), and the grey shaded zone is the total possible area we derive if we let all the parameters vary within their error boxes. The points corresponding to the fits of R96, MC02 and this work (Fig. 2, point B03) are superimposed to the plot.

We first note that the orbital elements of the orbit are rather badly constrained by the observations, as they are very sensitive to the measurement uncertainties. In fact, they are especially sensitive to the errors on ${\rm d}p/{\rm d}t$, ${\rm d}\psi /{\rm d}t$, and d. It also appears that giving independent error bars for aand e is impossible, as their variations are coupled: a can be as high as 45 AU but in this case e must fall around 0.1, while if a=27 AU, e must range between 0.4 and 0.55.

It could also appear surprising that the "B03'' point does not appear in the right plot in the middle of the black error box corresponding to d=140 pc. Indeed, this box is the exact counterpart in (a,e) space to the error box centered around the "B03'' point in $({\rm d}p/{\rm d}t,{\rm d}\psi/{\rm d}t)$ space in the left plot. This is in fact an illustration of the non-linearity of the relationship within the error box. If starting from the "B03'' point in the left plot, we let ${\rm d}p/{\rm d}t$ vary, we see that epasses through a minimum within the error box, and the starting point is in fact close to the minimum. The left edge of the box in the right plot corresponds in fact to that physical minimum, and not to the limits of the error box. This is the reason why the "B03'' point appears close to the left edge of the box in the right plot. In order to better emphasize this point, the edges of the outlined boxes in the right plot have been drawn in thick lines when they correspond to a physical limit, and in thin lines when they just correspond to the limits of the error box in $({\rm d}p/{\rm d}t,{\rm d}\psi/{\rm d}t)$ space.

3 A circumbinary disk dynamical study

As suggested above, we conjecture that the ring shape of the circumbinary disk around GG Tau A is a consequence of a dynamical sculpting by the inner orbital motion of GG Tau A, but also by tidal interaction with the outer binary GG Tau B. In fact, we roughly expect the inner edge of the disk to be fixed by the orbital motion of GG Tau A, and the outer edge to result from the interaction with GG Tau B. Now, before carrying out a dynamical study of the full system (Beust & Dutrey 2005), we first want to perform preliminary studies with a simplified description in order to better analyse the role of each parameter. In this study then, we restrict ourselves to the study of GG Tau A itself surrounded by its circumbinary disk. What we want to check here is whether the orbital solution derived above (Eq. (3)) is compatible with the observed inner edge of the disk at $R_{in} = 180\pm5~\mbox{AU}\times(d/140~\mbox{pc})$.

We are thus back to a classical circumbinary disk problem. Circumbinary disks are known to be cleared from inside by interaction with the components of the binary, leading to a sharp inner edge located well outside the orbit of the central binary. There is observational evidence for this fact, and this was predicted by many theoretical studies (Artymowicz & Lubow 1996,1994; Bate 2000). All these simulations show that the inner gap of the disk opens within a few orbital periods of the binary, with spiral density waves extending far in the remaining disk. Artymowicz & Lubow (1996) show that some of the tidally eroded mass can flow towards the individual stars by so-called "streamers''.

An interesting outcome of this theory is the ratio of the radius of the inner gap of the disk to the semi-major axis of the binary. This ratio depends on the eccentricity of the orbit, but Artymowicz & Lubow (1994) showed that it should fall between $\sim $2 and $\sim $4, and may be closer if a significant viscosity is present in the disk. These results were recently confirmed by symplectic integrations using the HJS software (Beust 2003). Another key parameter for this is the mass ratio $\mu= M_2/(M_1+M_2)$ of the secondary star (of mass M2) to the whole mass of the binary M=M1+M2. The mass ratio $\mu$ ranges from 0 if the secondary mass is negligible to 0.5 for two stars of equal masses. Artymowicz & Lubow (1994) and Beust (2003) present results for $\mu=0.1$ and $\mu=0.5$. Here we want to perform a specific study for GG Tau A. We thus make a similar study to that of Beust (2003), using the HJS symplectic software, and taking for $\mu$ a value corresponding to GG Tau A, i.e. $\mu=0.47$ if we assume the mass determinations of White et al. (1999) (see also Table 2). For this kind of study, SPH and symplectic integrations lead to similar results. Each method has its own advantages. SPH allows one to take gas dissipation into account but is rather limited in integration time, while symplectic methods allow very long-term integrations thanks to a long time-step, but cannot handle gas dynamics. As shown in Beust (2003), for a standard problem like the inner edge sculpting of a circumbinary disk, both methods lead to the same output, since the gravitational processes largely dominate the dynamics close to the edge of the disk. However, long term symplectic integrations presented in Beust (2003) reveal that after a few hundreds of orbital periods (i.e., the typical integration time for SPH runs), the sculpting process of the disk is not fully achieved. This motivates the use of our HJS integrator here.

Note that the problem is scale invariant, all distances rescaling with the semi-major axis a of the binary. We perform various integrations with HJS for different values of the eccentricity e. As in Beust (2003), the initial disk contains 105 particles located between r=1.5a and r=3.5a from the center of mass of the system. We take an inner value for r=1.5abecause due to the binarity, particles within one a to the stars are ejected (and lost to the simulation) in less than a few orbits. Following standard consideration on disks, we assume a surface density $\propto$r-1 for the particle distribution. The particles are randomly chosen with eccentricities between 0 and 0.1 and inclinations with respect to the binary orbital plane ranging between 0 and  $3\hbox {$^\circ $ }$.

The result is shown in Fig. 3. We note the similarity of this study to the preceding results. Depending on the eccentricity, the inner edge of the disk ranges between 2a and 3.3a. As the inner edge of the disk is known, this provides constraints on the orbit of the binary that we want to compare to those deduced from astrometry. More specifically, given an eccentricity value e, one can compute the semi-major axis range corresponding to an inner edge at the observed location. Letting e vary between 0 and 1 subsequently defines a band-shaped zone in (a,e) space compatible with the observed inner edge. It is then of interest to compare it to the error box we deduced from the fit of the orbital motion (Fig. 2).

This comparison is done in Fig. 4. The band we deduce from the location of the inner edge of the disk is superimposed on a plot similar to the right plot of Fig. 2. However, that band implicitly depends on the distance d via the measurement of the inner edge of the disk. It is then better to fix the value of the distance d to make the comparison. In Fig. 4, d is fixed to 140 pc, but similar plots for other values of d can be made.


  \begin{figure} \par\includegraphics[angle=-90,width=8.8cm]{2441f3} \end{figure} Figure 3: Approximate radius of the inner edge of the circumbinary disk as computed with HJS integrations, in units of the semi-major axis a of the binary, for various values of e, and $\mu =M_2/(M_1+M_2)=0.467$ corresponding to GG Tau A, taken from White et al. (1999).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm]{2441f4.eps} \end{figure} Figure 4: Comparison of the error boxes for the orbit of GG Tau A in (a,e) space deduced from astrometry (dark grey) and circumbinary disk dynamics (light grey). The conventions for this plot are the same as for the right plot of Fig. 2, except that the distance d has been fixed here to 140 pc for clarity. The box outlined in black is the error box for a and e we derive from the fit of the astrometric data if we let ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$ vary within their error bars ($1 \sigma $ coming from the astrometric data). The dark grey surrounding zone is the one we derive if we take into account the uncertainties on all the other parameters. The light grey band is the compatibility zone with the observed location of the inner edge of the disk. It is computed assuming that the location of the inner edge is known with $\pm $5 AU accuracy, but also taking a $\sim $$2\%$ uncertainty on the theoretical edge locations shown in Fig. 3.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm]{F5A.eps} \end{figure} Figure 5: Same as Fig. 4, but now the error bars on ${\rm d}p/{\rm d}t$, ${\rm d}\psi /{\rm d}t$, and on the location of the observed inner edge of the circumbinary disk have been increased by a factor of 3. There is now an intersection zone between the two boxes (outlined in white).
Open with DEXTER

The major outcome of Fig. 4 is that the two compatibility zones do not intersect. There is thus a discrepancy (already noted in MC02) between the astrometric fit of the orbital motion of the binary and the size of the disk gap. This discrepancy might not be real, and could be due to the difficulty to put relevant constraints on  ${\rm d}p/{\rm d}t$. As a matter of fact, Fig. 1 obviously shows that the error bars on the measurements of the separation are underestimated, otherwise it would not be compatible with any orbital motion.

We therefore decided to perform a new fit, but arbitrarily enlarging the error bars on ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$ by a factor of 3. This is equivalent to performing a $3 \sigma $ fit instead of $1 \sigma $ in Fig. 4. The results, shown in Fig. 5, follow the same plotting conventions as Fig. 4. Note that in Fig. 5 we also increase the error bar on the observed location of the inner edge of the circumbinary disk by a factor of 3, bringing it up to $\pm $15 AU. The difference to Fig. 4 is striking. There is now a zone, outlined in white in the figure, and centered around $a\simeq 62~$AU and $e\simeq 0.35$, that is compatible with all the constraints within the error bars.

Finally, in Fig. 6, we still assume those lower limits, but we also assume that the errors bars on the measurement of the inclination $\lambda $ and position angle $\phi $ of the circumbinary disk given by GDS99 may be underestimated; we thus multiply them by 3, fixing them to $\pm $ $3\hbox {$^\circ $ }$ and $\pm $ $6\hbox {$^\circ $ }$ respectively. In the left plot of Fig. 6, only the error bar on $\lambda $ is enlarged, and in the right plot only that on $\phi $. In a coplanar system, the position angle $\phi $ of the disk (as defined here) is also the projection of the rotation axis of the system onto the plane of the sky. The binary is currently viewed at position angle  $\psi\simeq 0$ with respect to the north (see appendix). Actually the real physical parameter is the difference of position angle between the disk and the binary.

Of course, when passing from Figs. 5 to 6, the error box in (a,e) space gets larger, but the effect is minor, and the main result noted in Fig. 6 holds: there is a more or less narrow intersection region centered around ( $a\simeq 62~$AU, $e\simeq 0.35$) that is compatible with all the observational constraints and with those placed by the disk dynamics. The size of the error box turns out to be essentially controlled by the accuracy of the astrometric measurements of the binary, much more than the other parameters.

Changing the error bars on the measurements of ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$ slightly changes the central fit values for these quantities. This subsequently changes the orbital fit to a and e. In Figs. 5-6 the central value is now close to a=42 AU and e=0.06. This is significantly different from Eq. (3) though still within the error box, and illustrates the extreme sensitivity of the orbital fit to the data.

  \begin{figure} \par\includegraphics[width=18cm]{F6A.eps} \end{figure} Figure 6: Same as Fig. 5 ($3 \sigma $ errors on astrometric data), but now the errors on the measurement of the inclination $\lambda $ ( left plot) and of the position angle $\phi $( right plot) of the circumbinary disk are enlarged by a factor of 3 with respect to the values of GDS99, setting them to $\pm $ $3\hbox {$^\circ $ }$ and $\pm $ $6\hbox {$^\circ $ }$ respectively. This shows that the astrometric errors are the dominant ones.
Open with DEXTER

Conversely, if we take a=62 AU and e=0.35, we can compute the corresponding rates of change ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$. We find

 \begin{displaymath}\frac{{\rm d}\psi}{{\rm d}t} = -1.423~\hbox{$^\circ$ }/\mbox{... ...\rm d}p}{{\rm d}t} = -4.36~\mbox{mas}~\mbox{yr}^{-1} \times d, \end{displaymath} (4)

which is still within the error bars. We note that the relative change with respect to Eq. (2) concerns mainly ${\rm d}p/{\rm d}t$, showing that this quantity is very poorly constrained. We also derive for that orbit choice $m=-12\hbox{$^\circ$ }$, showing that the binary is now closer to periastron than in the previous fit.

4 Discussion

The preceding study reveals a discrepancy between the astrometric fit of the orbit of GG Tau A and the size of the internal gap of its circumbinary disk. We have two possible choices for the orbit of GG Tau A: an orbit resulting from the fit of the astrometric data, characterized by (a=32.4 AU, e=0.34), hereafter referred to as Orbit A1, but that is incompatible with the location of the inner edge of the disk, and another orbit, twice as large, characterized by (a=62 AU, e=0.35), hereafter referred as Orbit A2, compatible with all constraints, but only marginally with the astrometric data.

Of course intermediate choices can be made, but our following discussion will focus on these two extreme cases. We may summarize the questions raised by the two possibilities for the dynamical status of GG Tau A and its disk:

We may now suggest the following ways to solve the discrepancy:
1.
The orbit is presently close to Orbit A1, but it used to be different (wider) in the past. In this case, obviously a formerly wider orbit should have cleared out the circumbinary disk, possibly up to 180 AU. Indeed, if we remember that GG Tau A is part of a quadruple system with GG Tau B, we may suggest that the gravitational interaction with GG Tau B could have caused a secular evolution of the orbit of GG Tau A. However, it is well known that in non-resonant situations, the semi-major axes are secular invariants, unless close encounters occur. Due to a smooth secular interaction with GG Tau B, we expect the eccentricity of GG Tau A to possibly evolve, but probably not its semi-major axis. In Fig. 4, we see that starting from the fit point B03 and letting only e evolve means moving along a vertical line, and that it cannot cross the gray band compatible with the disk gap this way. This conclusion needs to be confirmed by a dedicated dynamical study, but it will probably hold.
2.
Another possibility is that the clearing process of the inner edge of the disk, combined with the possibly still active accretion process onto the binary (via streamers), leads to a viscous interaction between the binary and the disk, resulting in a decrease of the semi-major axis of the orbit. According to this picture, the binary would transfer angular momentum to the inner edge of the disk and undergo a kind of inward orbital migration such as described for protoplanets by Ward (1997) and Trilling et al. (1998). What is not clear with this process is why the inner edge of the disk would not similarly migrate inward together with the orbit due to viscosity, keeping the ratio inner edge/semi-major axis constant, close to the nominal value we derive in this paper. The only way to solve this discrepancy is to suppose that the inner edge of the disk has been stabilized for a long time, but that accretion onto the binary still keeps going on thanks to spiral density waves launched across the disk by successive periastron passages of the outer pair GG Tau B, a process that clearly acts on a much longer timescale. In that case the orbit of GG Tau A could still be shrinking without affecting the inner edge of the disk. Of course this is very speculative and needs a dedicated study. Actually the clearing process of the inner disk may be not achieved yet, so that the equilibrium results presented in Fig. 3 could not apply. However, the calculations presented in Beust (2003) show that a circumbinary gap is never wider during the clearing phase than at the final equilibrium.

3.
The orbit is close to Orbit A1, but the size of the inner gap is not due to the sole interaction with the binary, but rather to another unseen (planetary-sized?) body orbiting the binary around $\sim $140 AU. This hypothesis, purely speculative, raises a few interesting questions. The minimum mass required for such a (massive) body to efficiently clear out the disk from $\sim $90 AU up to $\sim $180 AU can be estimated from previous work on gap opening processes by planets in protoplanetary disks. The approximate width of the gap cleared by a protoplanet in a disk is theoretically given by Takeuchi et al. (1996):

 \begin{displaymath}\Delta r \simeq 0.9 a\left(\frac{q^2}{\alpha h^2}\right)^{1/3}, \end{displaymath} (5)

where q is the ratio of the planet mass to the central mass, his the dimensionless scale height of the disk, $\alpha$ is the dimensionless viscosity parameter, a is the orbital semi-major axis of the planet. The validity of this estimate has been confirmed by various numerical studies (Nelson & Benz 2003; Trilling et al. 1998). Following Takeuchi et al. (1996), we take h=0.15 (from GDS99) and $\alpha=10^{-2}$ (as a standard value for such disks). We assume a=135 AU in order to put the planet in the middle of the region to be cleared. We want to have $\Delta r\simeq 90~$AU in order to clear the disk between 90 AU and 180 AU. With the above quoted numerical values, this requires a body of mass $\sim $12.8 Jupiter masses ( $0.012~M_\odot$). Equation (5) is theoretically valid for small gaps; the simulations by Trilling et al. (1998) show that for larger gaps such as required here, the gap opened by a given planet is sightly smaller than predicted by Eq. (5). Hence the required planetary masses may be slightly larger than the above estimates. From Fig. 9 of (R96) which gives the J magnitude inside the hole of the GG Tau ring and using the evolutionnary tracks from (Baraffe et al. 1998), we cannot rule out the presence of an unseen object up to a mass of $0.02~M_\odot$. Note however, that the orbital stability of such a body needs also to be addressed.
4.
The orbit is actually close to Orbit A2, but the disk and the orbit of GG Tau A are not coplanar. Up to now, we have always assumed that the midplane of the disk coincides with the orbital plane of GG Tau A, and this was used by all authors to constrain the orbit of the binary (R96; MC02; Tamazian et al. 2002). Dropping the coplanar assumption means in our formalism releasing all constraints on the angles $\phi $ and $\lambda $. More specifically, $\phi $ and $\lambda $ are defined relative to the image of the circumbinary disk and remain unchanged. We shall call now $\lambda'$ the inclination with respect to the plane of the sky of the orbital plane of the binary and $\phi'$ the position angle (with respect to north) of the projection of the rotation axis of its orbit onto the plane of the sky. Of course if the disk and the orbit are coplanar we have $\lambda=\lambda'$ and $\phi=\phi'$, but this is not true in general. If we assume for a and e the values corresponding to Orbit A2, and if we assume for a and ethe fit (2), it is possible to solve the geometrical and kinematic Eqs. (A.7)-(A.11) (see Appendix) for $\lambda'$ and $\phi'$, i.e., finding an orbital plane orientation compatible with all the constraints. Performing this, we find one possible value for the orbital inclination $\lambda'$ with respect to the plane of the sky, namely

\begin{displaymath}\lambda'=125.4\hbox{$^\circ$ }, \end{displaymath} (6)

and 4 possible values for $\phi'$:

\begin{displaymath}\phi'\!=\!24.3\hbox{$^\circ$ },\:\phi'\!=\!-155.7\hbox{$^\cir... ...hbox{$^\circ$ },\: \mbox{and}~\phi'\!=\!164.4\hbox{$^\circ$ }. \end{displaymath} (7)

For each set of values ( $\lambda',\phi'$), it is then possible to compute the tilt angle $\tau$ between the orbital plane and the midplane of the disk. We find

\begin{displaymath}\tau =90.3\hbox{$^\circ$ },\; \tau=21.4\hbox{$^\circ$ },\; \tau=23.7\hbox{$^\circ$ }\;\mbox{and}\; \tau=89.4\hbox{$^\circ$ }. \end{displaymath} (8)

The disk appears then either inclined by $\sim $ $20\hbox{$^\circ$ }$ with respect to the orbital plane of GG Tau A, or nearly perpendicular to it. A non-coplanar circumbinary disk could appear surprising, though it should not be immediately ruled out. Disks are sometimes non-coplanar. Koresko (1998) reports for example the detection of a possibly non-coplanar circumstellar (not circumbinary) disk orbiting HK Tauri B. From a theoretical point of view, Erwin (1998) shows that some orbits in circumbinary disks can be destabilized in the direction of offplane motion. However, in such a case we would expect the disk to achieve a complex, highly warped structure, which does not seem to be the case here. If the disk is actually non-coplanar, this should be probably related to tidal interactions with the outer binary GG Tau B. This issue is dynamically investigated in Beust & Dutrey (2005), and the result confirms what is suggested above. Assuming non-coplanarity leads to a disk that is much thicker than observed.
5.
Finally, the orbit is close to Orbit A2, the disks are almost coplanar but the error bars on the astrometry deduced from optical/NIR measurements are strongly underestimated.

5 Conclusion

Using the astrometric data of the binary GG Tau A available in the literature, we derive an accurate fit of its orbit (a=32.4 AU, e=0.34), consistent with recent similar observational works. However, the inner radius of the circumbinary ring surrounding GG Tau A, measured from both NIR and millimetric observations is dynamically incompatible with the fit of the orbit provided by the astrometric data. Our dynamical study shows that the ring inner gap is approximately twice as large as should be expected if we assume the orbit given by the astrometric fit. To explain the discrepancy, we introduce larger error bars i) on the astrometric data and ii) on the disk PA and inclination, these effects being separately investigated. Larger error bars on the disk PA and inclination do not significantly affect the fit of the orbit. By introducing larger astrometric error bars, we find another possible orbit (a=62 AU, e=0.35) which is fully compatible with the observed inner gap, but only marginally consistent (3$\sigma $) with the astrometric data. There are four possible ways to solve this discrepancy: i) the orbit of GG Tau A has undergone a significant secular evolution in the past, due to the interaction with GG Tau B; ii) a large (and massive) circumbinary planet orbits GG Tau A around $\sim $140 AU clearing the disk further out than normally occurs; iii) the disk and the binary are not coplanar; or iv) more likely the system is essentially coplanar but the error bars on the astrometric data are underestimated.

All these hypotheses must be dynamically investigated, in a framework where the dynamics of the quadruple stellar GG Tau system and disk should be treated as a whole. This is the subject of a forthcoming paper (Beust & Dutrey 2005).

Acknowledgements
We thank Gaspard Duchêne, Stéphane Guilloteau and Mike Simon for many fruitful discussions about the GG Tau system. All the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI). Figures were done using the GILDAS software.

Appendix A: Deriving orbital elements from astrometric measurements of the binary

The orbital plane of the binary is assumed to be inclined by angle $\lambda $ with respect to the plane of the sky. Let us suppose that the radius vector of the binary achieves the polar coordinate set $(r,\theta)$ in its plane; r is the physical separation between the two stars and $\theta$ is a polar angle, $\theta=0$ corresponding to the radius vector aligned with the apparent small axis of the image of the disk. This small axis is tilted by position angle $\phi $ with respect from north. The components (px,py) of the projected radius vector, once expressed in an ( $\vec{I},\vec{J}$) referential frame of the plane of the sky, where $\vec{I}$ points towards west and $\vec{J}$towards north (see Fig. A.1), read

\begin{displaymath}\left\{\begin{array}{rcl} p_x & = & r\left(\cos\theta\cos\lam... ...s\lambda\cos\phi+\sin\theta\sin\phi\right) \end{array}\right.. \end{displaymath} (A.1)

We immediately derive the projected separation p:

\begin{displaymath}p=\sqrt{p_x^2+p_y^2}=r\sqrt{1-\cos^2\theta\sin^2\lambda}. \end{displaymath} (A.2)

Similarly, we get the components of the apparent velocity

\begin{displaymath}\left\{\begin{array}{rcl} v_x & = & \displaystyle \frac{\part... ...n\theta\cos\phi\cos\lambda\right)v_\theta \end{array}\right. , \end{displaymath} (A.3)

where vr and $v_\theta$ are the radial and orthoradial velocities respectively in the orbital motion.

We also derive the rate of change of the projected separation pand of the position angle $\psi $ of the radius vector

                                          v$\displaystyle \frac{{\rm d}p}{{\rm d}t}$ = $\displaystyle \frac{p}{r}~v_r+\frac{\partial p}{\partial\theta} \frac{v_\theta}... ...os\theta\sin\theta\sin^2\lambda}{\sqrt{1-\cos^2\theta\sin^2\lambda}} ~v_\theta;$ (A.4)
$\displaystyle \frac{{\rm d}\psi}{{\rm d}t}$ = $\displaystyle \frac{p_xv_y-p_yv_x}{p^2}=\frac{\cos\lambda} {1-\cos^2\theta\sin^2\lambda}~\frac{v_\theta}{r}\cdot$ (A.5)

Now, the binary is currently seen at position angle  $\psi\simeq 0$. Demanding $\psi=0$ means setting px=0. This fixes the value of $\theta$. Once this is done is then possible to replace  $\sin\theta$ and  $\cos\theta$ by their corresponding values in the expressions of p, ${\rm d}p/{\rm d}t$ and ${\rm d}\psi /{\rm d}t$. After some straightforward algebra, we get
  
                                     p = $\displaystyle r~\frac{\vert\cos\lambda\vert}{\sqrt{1-\sin^2\lambda\sin^2\phi}} ;$ (A.6)
$\displaystyle \frac{{\rm d}\psi}{{\rm d}t}$ = $\displaystyle \frac{1-\sin^2\lambda\sin^2\phi}{\cos\lambda} ~\frac{v_\theta}{r} ;$ (A.7)
$\displaystyle \frac{{\rm d}p}{{\rm d}t}$ = $\displaystyle \frac{\vert\cos\lambda\vert} {\sqrt{1-\sin^2\lambda\sin^2\phi}}~v... ... sign}(\cos\lambda)\sin^2\lambda} {\sqrt{1-\sin^2\lambda\sin^2\phi}}~v_\theta ,$ (A.8)

which provides the relationship between the observed quantities ( $p, {\rm d}p/{\rm d}t, {\rm d}\psi/{\rm d}t$) and the local kinematic quantities ( $r,v_r,v_\theta$).


  \begin{figure} {\includegraphics[angle=-90,width=0.7\columnwidth]{2441fA1} }\end{figure} Figure A.1: Projection of the circumbinary disk (grey ring) and of the GG Tau A binary (black spheres) onto the plane of the sky, showing the definition of the vectors ( $\vec{I},\vec{J}$), the position angles $\psi $ and $\phi $, and the projected separation p. The disk itself is inclined by oriented angle $\lambda $ with respect to the plane of the sky.

As the observation shows obviously that ${\rm d}\psi/{\rm d}t<0$, this shows that $\cos\lambda<0$ ( $v_\theta>0$ by definition), meaning as noted in GDS99 that the disk is viewed from the south side (actually $\lambda=180\hbox{$^\circ$ }{-}37\hbox{$^\circ$ }$).

Now, r, vr and $v_\theta$ are related to the orbital elements by a classical Keplerian formalism, namely:

 
                      r = $\displaystyle \frac{a(1-e^2)}{1+e\cos v} ;$ (A.9)
vr = $\displaystyle \frac{e\sin v}{\sqrt{1-e^2}}~\sqrt{\frac{GM}{a}} ;$ (A.10)
$\displaystyle v_\theta$ = $\displaystyle \frac{1+e\cos v}{\sqrt{1-e^2}}~\sqrt{\frac{GM}{a}} ,$ (A.11)

where M is the total mass of the binary, a is the semi-major axis of the orbit, e is its eccentricity, and v is the current true anomaly. It is then only a matter of algebra to combine the two preceding sets of equation and finally express the orbital elements as a function of the observed quantities. The semi-major axis a finally reads (assuming $\cos\lambda<0$)
                                  $\displaystyle \frac{1}{a}$  =  $\displaystyle \frac{\sin^2\phi^2\sin^2\lambda-1}{GM\cos^2\lambda} \left(\frac{{... ...lambda}{GM\cos^2\lambda} p\frac{{\rm d}\psi}{{\rm d}t}\frac{{\rm d}p}{{\rm d}t}$  
    $\displaystyle -\frac{1-\cos^2\phi\sin^2\lambda}{GM\cos^2\lambda} p^2\left(\frac... ...ght)^2\!-\frac{\cos\lambda}{\sqrt{1-\sin^2\phi\sin^2\lambda}}\frac{2}{p} \cdot~$ (A.12)

The eccentricity is less easy to express, but it is related to the semi-major axis by the quantity a(1-e2), which achieves a very simple form:

\begin{displaymath}a(1-e^2)=\frac{p^4}{GM\cos^2\lambda}~\left(\frac{{\rm d}\psi}{{\rm d} t}\right)^2 \cdot \end{displaymath} (A.13)

This quantity is actually related to the specific angular momentum $C=\sqrt{aGM(1-e^2)}$; interestingly, it does not contain  ${\rm d}p/{\rm d}t$. As ${\rm d}p/{\rm d}t$ is by far the most weakly constrained observational quantity (see text), we therefore have here a relation between a and e that we hope to be better constrained than a and e themselves.

We finally derive the true anomaly v as:

\begin{displaymath}\left\{\begin{array}{rcl} e\cos v & = & \displaystyle -1\!-\!... ...si}{{\rm d}t}\frac{{\rm d}p}{{\rm d}t} \end{array}\right.\cdot \end{displaymath} (A.14)

The true anomaly is then related to the mean anomaly m by standard Keplerian formalism.

References

 
Copyright ESO 2005