# Dynamics of the young multiple system GG Tauri II. Relation between the stellar system and the circumbinary disk

H. Beust, Anne Dutrey

► To cite this version:

H. Beust, Anne Dutrey. Dynamics of the young multiple system GG Tauri II. Relation between the stellar system and the circumbinary disk. Astronomy & Astrophysics - A&A, 2006, 446, pp.137-154. 10.1051/0004-6361:20053163. [hal-00023656v1](https://hal.science/hal-00023656v1)

HAL Id: hal-00023656

https://hal.science/hal-00023656v1

Submitted on 13 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Dynamics of the young multiple system GG Tauri

## II. Relation between the stellar system and the circumbinary disk*

H. Beust<fnref n="1" /> and A. Dutrey<fnref n="1" />,<fnref n="2" />

Received 31 March 2005 / Accepted 15 August 2005

## ABSTRACT

The quadruple young system GG Tauri is an example of a multiple T Tauri system. Its consists of two binaries, the brighter one (GG Tau A) being surrounded by a ring-shape circumbinary disk. In a recent paper, we performed a dynamical study of circumbinary and showed that there is an apparent discrepancy between the orbital fit of GG Tau A and the observed inner edge of the disk. In this paper, we now investigate the dynamics of the whole quadruple system together with the disk. We show that it is possible to design an orbital configuration between the two binaries in such a way that the outer profile of the circumbinary ring may be explained by tidal interaction with repeated periastron passages of the outer binary GG Tau B. We show that the observed characteristics of the disk are not compatible with some orbital configurations, such as those giving rise to the Kozai resonance. Surprisingly, the outer binary GG Tau B appears only marginally stable against tidal disruption by GG Tau A. GG Tau B appears stable only if its internal orbit is retrograde with respect to the motion of its center of mass around GG Tau A. We also find that the CB disk should be almost coplanar with the inner binary.

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

## 1. Introduction

GG Tauri is a quadruple T Tauri systems consisting of two binaries. The first pair (GG Tau A) is brighter and closer (0.25"). The second pair (GG Tau B) is wider (1.48") and located 10.1" to the south (Guilloteau et al. 1999, hereafter GDS99; White et al. 1999). GG Tau A is well known to be surrounded by a circumbinary dusty and gaseous disk that was spatially resolved in both millimeter and near infrared wavelength domains (Dutrey et al. 1994; GDS99). This disk is the prototype example of a circumbinary disk in a young stellar system.

This disk has been the subject of many investigations the past. The disk is ring-like (GDS99) with well defined edges. The continuum dust data reveal that almost 70% of the material is located between 180 and 260 AU around GG Tau A, while the ¹³CO(2–1) line emission shows that the rest of the material extends up to 800 AU or more. The total mass of the circumbinary material (H₂ + dust) is ~0.12 M☉ (GDS99).

The Keplerian motion of this disk around the GG Tau A pair has been unambiguously identified. The first indication of a velocity gradient in the disk was given by Kawabe et al. (1993),

furthermore GDS99 made a careful analysis of the velocity map of the disk in ¹²CO(2–1) and ¹³CO(2–1) emissions and conclude that the motion of the disk was essentially Keplerian around a central mass of 1.28 ± 0.07 M☉. Mass determinations for all four components of the system are given by White et al. (1999), derived from stellar evolution models. They give masses of 0.78 ± 0.1 M☉ and 0.68 ± 0.03 M☉ for GG Tau Aa and Ab, and 0.12 ± 0.02 M☉ and 0.044 ± 0.006 M☉ for GG Tau Ba and Bb. Note that the total mass inferred for GG Tau A remains compatible with the dynamical determination of GDS99 within the error bars, and that GG Tau B is much less massive than GG Tau A. GG Tau Bb seems even to be a substellar object. For our simulations, we adopt these mass determinations.

The relative motion of the components of GG Tau A has been now monitored for more than 10 years (see Table 1 from Beust & Dutrey 2005, hereafter BD05), providing additional kinematic constraints. Combined with the assumption that the circumbinary disk is coplanar with GG Tau A, this allows one to derive a full orbital solution for GG Tau A. This was done by Roddier et al. (1996); McCabe et al. (2002); Tamazian et al. (2002). In a recent paper (BD05), we reinvestigate this issue using the most recent astrometric data, and derive in this way

<footnote n="1">¹ Laboratoire d'Astrophysique de Grenoble, Université J. Fourier, BP 53, 38041 Grenoble Cedex 9, France
e-mail: herve.beust@obs.ujf-grenoble.fr</footnote>

<footnote n="2">² Laboratoire d'Astrodynamique, d'Astrophysique et d'Aéronomie de Bordeaux, 2 rue de l'Observatoire, BP 89, 33270 Floirac, France</footnote>

<footnote>* Appendix A is only available in electronic form at
http://www.edpsciences.org</footnote>
an orbit characterized by a semi-major axis *a* = 32.4 AU and an eccentricity *e* = 0.34. This is in agreement with other recent work (McCabe et al. 2002; Tamazian et al. 2002). In BD05, we also perform dynamical simulations of the circumbinary disk on the basis of this orbital solution, and show that this orbit is not compatible with the location of the inner edge of the circumbinary disk observed at 180 AU, if we suppose that the truncation results from the sole interaction with the binary; it should be expected to be closer to the star by a factor of 2. The observational data about this inner edge location are robust up to ~10% and come from independent measurements obtained at several wavelengths (NIR, millimeter) (Dutrey et al. 1994; GDS99; Roddier et al. 1996). We can however solve this discrepancy if we assume that the error from bars on the astrometric data have been severely underestimated.

The aim of paper is to further investigate the dynamics of
this system, but considering now the whole quadruple system.
We adopt a didactic approach, where we describe the physics
related to the dynamics of multiple systems in general. These
papers are a first attempt to provide a link between various
observations of a given astrophysical object at several wave-
lengths, and dynamics. This is why we first describe 3-body
dynamics where GG Tau B is considered as a single object, be-
fore reintroducing its binary nature.

We first discuss the general shape of the disk. While it is
reasonable to think that the inner edge of the disk is due to tidal
interaction with the inner binary GG Tau A (despite the above
noted discrepancy), the outer profile needs some other expla-
nation. As noted by GDS99, there is no way to fit the outer
profile with a single power law. It rather looks like a sharp
ring ending at 260 ± 5 AU superimposed onto a more tenu-
ous tail following a power law of index s = -2.75. It is thus
tempting to attribute this sculpting of the outer edge to an in-
teraction with the outer binary GG Tau B. We investigate this
issue in the present paper. Our second goal is to investigate the
secular evolution of the quadruple system, as past secular evo-
lution could help to solve the discrepancy between the orbit of
GG Tau A and the inner edge of the disk. In Sect. 2, we perform
three-body simulations, treating GG Tau B as a single body. We
show that it is possible to explain the general shape of the cir-
cumbinary disk by tidal interaction with GG Tau B, provided
the wide orbit of GG Tau B around GG Tau A has some char-
acteristics we list. Then in Sect. 3, we reintroduce the binary
nature of GG Tau B. We show that this has only little effect on
the sculpting of the disk, but we also show that the GG Tau B
pair is only marginally stable against tidal disrupting due to
GG Tau A. GG Tau B appears stable only if the relative mo-
tion between its components is retrograde with respect to the
wide orbit between the two pairs. In Sect. 4, we discuss the
relevance of our calculations, and show that the discrepancy
between the inner edge of the disk and the orbital solution for
GG Tau A remains unsolved. We give possible clues for solv-
ing this problem, and we conclude in Sect. 5. by deriving orbital
solutions (Roddier et al. 1996; McCabe et al. 2002; Tamazian
et al. 2002).

The dynamical simulations presented in this paper, as well
as those of BD05, are performed with the HJS symplectic in-
tegrator (Beust 2003). HJS (Hierarchical Jacobi Symplectic;

Beust 2003) is a variant of the popular symplectic integra-
tion method WHM (Wisdom-Holman Mapping; Wisdom &
Holman 1991; Levison & Duncan 1994; Duncan et al. 1998),
but designed for the dynamics of hierarchical stellar *N*-body
systems, while the original method accounts for planetary sys-
tem dynamics. Symplectic integration allows us to integrate
very fast (with a large time-step) *N*-body systems for which
a dominant Keplerian part may be identified. This is the case in
planetary system dynamics but also for hierarchical systems.
Moreover, symplectic integrations schemes ensure that the en-
ergy error is bound (Chambers 1999).

**2. The three-body numerical study**

**2.1. General features**

In this section, we investigate the dynamics of the GG Tau A+B
system with the HJS software, but merging the two components
of GG Tau B into a single one (assuming its total mass is the
sum of the mass of each component, see White et al. 1999),
thus reducing the system to triple. Our aim is here to investigate
the perturbing action of GG Tau B on the circumbinary disk of
GG Tau A, which is mainly due to the total mass of GG Tau B
rather than its internal structure. The binary nature of GG Tau B
will be reintroduced later.

Performing a run under this assumption means choosing
two initial orbital configurations for the system: the mutual or-
bit of the two components of GG Tau A and the wider orbit of
the outer body GG Tau B around the center of mass of the cen-
tral pair. In the following, we refer to these orbits as “AA” and
“BA” respectively. Due to mutual gravitational perturbations,
these orbits are supposed to secularly evolve along any run we
perform. This evolution will be computed by the HJS integra-
tor. Hence the orbits we give below, which are designed to be
compatible with the *present* orbital configuration of GG Tau,
will be given as *initial* conditions to the runs. As the present
age of GG Tau is 1–2 Myr, it could appear more appropriate
to try to fit these conditions at that age after the beginning of
each run. We nevertheless think that this refinement is unnec-
essary, for several reasons. First it is very difficult to guess the
past dynamical evolution of GG Tau, as it is well known that
very young multiple systems are subject to violent dynamical
relaxation such as component ejection. Giving the present con-
figuration as the initial condition is not worse than integrating
backwards 1 or 2 Myr in the past to derive hypothetical ini-
tial conditions. Second, it is well known (unless when mean-
motion resonances are present, which is not the case here) that
phase angles like mean anomalies have a negligible role in the
secular evolution. Starting with the present values is equivalent
to starting with any other set of values. Third, even when Kozai
resonance is present (see below), the secular evolution of the
stellar system is moderate over a time span of 1 or 2 Myr. As
we are interested in mid-term stability criteria for the disk, the
result is expected to be roughly the same irrespective of the
starting point on the secular evolution.

The initial circumbinary disk itself will be simulated with a
set of non-interacting massless test particles, initially randomly
chosen, orbiting GG Tau A. The orbital eccentricities of the
particles are chosen between 0 and 0.1, and their inclinations with respect to their midplane between 0 and 3°. Depending on the runs, the initial disk holds 50 000 or 60 000 particles. Note that we choose to neglect the self gravity of the disk since the mass ratio between the disk and the total mass of the GG Tau A binary is ~0.1. Pierens et al. (2005), taking into account the self gravity of the GG Tau disk, reveal that its contribution is a second order effect mostly related to the mass inflow towards the central binary and likely does not affect the general dynamics of the system.

### 2.1.1. Orbit AA

Orbit AA is well constrained by the astrometric data. If we assume that the circumbinary disk of GG Tau A lies in the same plane as its orbit, then an astrometric fit of the projected motion of its two components is enough to fully characterize the orbit. In BD05, we derive in this way an orbit characterized by a semi-major axis $a = 32.4$ AU and an eccentricity $e = 0.34$. This is in agreement with other recent work (McCabe et al. 2002; Tamazian et al. 2002). However, we show in BD05 that this orbit is not compatible with the location of the inner edge of the circumbinary disk observed at 180 AU, if we suppose that the truncation results from the sole interaction with the binary; with $a = 32.4$ AU and $e = 0.34$ it should be expected two times closer to the star. This orbit will be referred as “Orbit AA1” in the following.

Conversely, if we let the errors bars on the astrometric data be larger than previously estimated, another orbit characterized by ($a = 62$ AU, $e = 0.35$) can be found, fully compatible with the observed location of the outer edge of the disk, while still within the enlarged error box deduced from the astrometric fit. This will constitute our “Orbit AA2”.

With Orbit AA2, the agreement between the shape of the disk and the astrometric data is still marginal. This is why in BD05, we drop the assumption that the circumbinary disk is coplanar with GG Tau A. It is then possible to constrain the orbital plane of GG Tau A in such a way that it is fully compatible with the astrometric data and with the inner edge of the disk at ~180 AU. 4 orbital solutions are found, hereafter referred to as “Orbits AA3–6”, two of them with a circumbinary disk inclined by ~20° with respect to Orbit AA, and nearly perpendicular to it for the two other solutions. The different orbits that we will consider for GG Tau A are summarized in Table 1.

### 2.1.2. Orbit BA

The relative orbit between the two binaries (BA) needs also to be specified. This orbit is much less well known than Orbit AA. All we know are the projected distances and position angles. The measured projected distance between the centers of mass of the two binaries is 10.4" (i.e. 1440 AU with a distance of GG Tau of 140 pc) towards the south from GG Tau A (PA = 180°). We have no kinematic data. Most of the parameters of that orbit will be free parameters in our study. If we take into account the observational constraints, only 4 free parameters are needed to fully characterize the orbit. These parameters will be

<table><caption>**Table 1.** The various initial configurations taken for the orbit of the GG Tau A binary (Orbit AA). The inclinations are given relative to the midplane of the circumbinary disk, and the mean anomaly corresponds to the present observed configuration of the binary.</caption><thead><tr><th>Orbit designation</th><th>Semi-major axis <i>a</i></th><th>Eccentricity <i>e</i></th><th>Inclination <i>i</i></th><th>Current mean anomaly <i>M</i></th></tr></thead><tbody><tr><td>AA1</td><td>32.4 AU</td><td>0.34</td><td>0°</td><td>179.9°</td></tr><tr><td>AA2</td><td>62.0 AU</td><td>0.35</td><td>0°</td><td>-20.8°</td></tr><tr><td>AA3</td><td>62.0 AU</td><td>0.35</td><td>90.3°</td><td>-56.0°</td></tr><tr><td>AA4</td><td>62.0 AU</td><td>0.35</td><td>21.4°</td><td>-56.0°</td></tr><tr><td>AA5</td><td>62.0 AU</td><td>0.35</td><td>23.7°</td><td>61.0°</td></tr><tr><td>AA6</td><td>62.0 AU</td><td>0.35</td><td>89.4°</td><td>61.0°</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** The radial profile of the initial disk common to all our simulation runs using Orbit A2 for GG Tau A: the surface density as a function of the distance *R* to the center of the disk.</figcaption></figure>

its periastron distance $q'$, its initial inclination $i'$ with respect to the midplane of the circumbinary disk, and two parameters relative to the present location of GG Tau B on its orbit: the mean anomaly $M'$ and the current altitude $z'$ above the plane of the disk along the line of sight (in the following, primed quantities will refer to Orbit BA). Alternatively, we may replace one of these parameters by the initial relative inclination $i_r$ between the two orbits, but the distinction holds only for situations where the disk and the central binary are not coplanar, otherwise $i'$ and $i_r$ coincide. Of course, not all combinations of these variables are possible.

As we will see below $q'$ and $i'$ (alternatively $i_r$) are by far the most important parameters controlling the overall dynamics of the system.

For runs using Orbit AA1 for GG Tau A, the disk initially holds 60 000 particles, with semi-major axes ranging from 70 to 1200 AU, chosen in such a way that it yields a surface density $\propto r^{-1}$. For the other runs, the disk holds 50 000 particles, with semi-major axes between 150 and 1200 AU, with the same surface density power index. The radial profile of that initial disk is shown in Fig. 1. Note that in both cases, the radial extension of the initial disk is significantly larger than for the real observed disk in both directions.
Starting from various initial orbital situations, we integrate the dynamics of the whole system with the HJS method over 15 Myr, assuming a time-step of 20 yrs (~1/20 of the smaller orbital period). In Beust (2003), it was shown that in general triple systems, assuming such a time-step in the HJS integrator, ensures stability of the energy of the whole system up to a few $10^{-5}$ in relative value. This is enough for our purpose. In the 4-body study (see below) we will need to reduce the time-step because some orbits are close to instability.

Note that the integration time is significantly larger than the estimated age of the GG Tau system. But our goal is to investigate the dynamical stability of the circumbinary disk and of the stellar system, so that integrating over a longer time-span is necessary.

For each run, we compare the characteristics of the observed disk to the ones we derive after gravitational sculpting by interaction with the stellar system. The most obvious criterion is the radial confinement. We must obtain a ring-like structure confined between ~180 AU and ~260 AU with a somewhat less sharp outer edge, and possibly an outer, more tenuous tail extending up to ~800 AU (GDS99). We may draw radial profiles for our simulated disks in order to compare them to the data. Note that in deriving these radial profiles, in order to achieve better resolution, we average all the remaining particles over their respective orbits: for each particle with its osculating orbit, we take into account 200 fictitious particles spread over the orbit following Kepler's second law. This applies to Fig. 1 and to all the profiles presented below. In each case we checked that if we built the profile only with the initial particles, we get the same, but more noisy profile. We use the information we have concerning not only the position of the particles, but also their orbits that are natural outputs of HJS computations.

Another output to be compared to the data is the vertical profile of the simulated disk. GDS99 give the scale height law of the disk:

$$H(r) = H_{300} \left( \frac{r}{300 \text{ AU}} \right)^h, \qquad (1)$$

with $H_{300} = 55 \pm 5$ AU and $h = 1.05 \pm 0.05$. Following again the rule that we want to use the information about the orbits of the particles, we adopt the following strategy: for each output, we compute the midplane of the disk of particles (by computing its mean angular momentum) and then we derive the mean orbital inclination $i_d$ of the particles with respect to that plane. At each distance we expect to have

$$H(r) \approx r \tan i_d. \qquad (2)$$

As in Eq. (1), the index $h$ is close to 1, we expect $i_d$ to be almost not dependent on $r$. With the quoted value of $H_{300}$, we expect $i_d \approx 10^\circ$.

## 2.2. Coplanarity of the disk with GG Tau A

In this section, we investigate the orbital configurations with a circumbinary disk initially coplanar with the central binary, i.e., runs with Orbits AA1 and AA2 from Table 1. We then test various orbital configurations for Orbit BA. Their list is given in Table 2.

<table><caption>**Table 2.** From BD05, the various initial configurations taken for Orbit BA in runs associated with a circumbinary disk coplanar with Orbit AA (Orbits AA1 and AA2) The initial inclinations are given relative to the midplane of the circumbinary disk, which in all these cases coincide with the relative inclination $i_r$ with respect to Orbit AA.</caption><thead><tr><th>Orbit designation</th><th>Periastron q'</th><th>Eccentricity e'</th><th>Inclination i'</th><th>Altitude z'</th><th>Current mean anomaly M'</th></tr></thead><tbody><tr><td>BA1</td><td>800 AU</td><td>0.5152</td><td>20°</td><td>2000 AU</td><td>150°</td></tr><tr><td>BA2</td><td>600 AU</td><td>0.4464</td><td>30°</td><td>600 AU</td><td>-150°</td></tr><tr><td>BA3</td><td>700 AU</td><td>0.3825</td><td>20°</td><td>600 AU</td><td>-150°</td></tr><tr><td>BA4</td><td>800 AU</td><td>0.3237</td><td>80°</td><td>600 AU</td><td>150°</td></tr><tr><td>BA5</td><td>1100 AU</td><td>0.1725</td><td>20°</td><td>600 AU</td><td>150°</td></tr><tr><td>BA6</td><td>800 AU</td><td>0.5152</td><td>160°</td><td>2000 AU</td><td>150°</td></tr><tr><td>BA7</td><td>600 AU</td><td>0.4464</td><td>160°</td><td>600 AU</td><td>-150°</td></tr></tbody></table>

<table><caption>**Table 3.** List of three-body runs with initial coplanar disk described in this paper. For each run, we list the initial orbital combination taken from Tables 1 and 2, the approximate resulting boundaries of the main surviving circumbinary ring, the tilt angle between the midplane of the disk and the orbital plane of the central binary and the mean inclination $i_d$ of the particles with respect to that midplane. We add a comment to make a first classification of the runs.</caption><thead><tr><th>Run #</th><th>Orbit combination</th><th>Ring boundaries (AU)<br>t = 2 Myr t = 15 Myr</th><th>Tilt angle i_d</th><th>Comment</th></tr></thead><tbody><tr><td>3b1</td><td>AA2+BA1</td><td>180-370 180-340</td><td>≤1° 6°</td><td>Reference</td></tr><tr><td>3b2</td><td>AA2+BA2</td><td>180-260 200-240</td><td>≤5° 20°</td><td>Reference</td></tr><tr><td>3b3</td><td>AA2+BA3</td><td>180-300 185-300</td><td>≤2° 7°</td><td>Reference</td></tr><tr><td>3b4</td><td>AA1+BA1</td><td>90-400 90-340</td><td>~1° 8°</td><td>Reference</td></tr><tr><td>3b5</td><td>AA2+BA4</td><td>190-360 205-360</td><td>~1° 20°</td><td>Kozai</td></tr><tr><td>3b6</td><td>AA1+BA4</td><td>95-300 110-300</td><td>~5° 20°</td><td>Kozai</td></tr><tr><td>3b7</td><td>AA2+BA5</td><td>180-550 180-530</td><td>~3° 10°</td><td>Large q'</td></tr><tr><td>3b8</td><td>AA2+BA6</td><td>180-500 180-500</td><td>~3° 20°</td><td>Retrograde</td></tr><tr><td>3b9</td><td>AA2+BA7</td><td>180-350 180-350</td><td>≤2° 20°</td><td>Retrograde</td></tr></tbody></table>

As suggested above the most important parameters are $q'$ and $i'$. We performed many runs that helped us to classify them. In this section, we only present a representative subset. The main characteristics and basic results of these runs are summarized in Table 3. All these runs lead to sculpting a ring-like circumbinary disk that more or less matches the observed one, and that remains stable for the time-span of the simulation (and probably indefinitely in most cases). This can be seen in the radial boundaries of the ring-like structures, shown at $t = 2$ Myr (the age of the system) after the beginning of the run, at $t = 15$ Myr (the end of the simulation). We also note that in all cases, the tilt angle between the disk and Orbit AA is very small, so that the disk remains coplanar with the central binary. The inclination dispersion $i_d$ must be compared to the expected ~10° deduced from the scale height of the disk. Obviously some runs fulfill this constraint while some others do not.

We associate each run with a comment such as "Reference", "Kozai"... (Table 3). Runs with the same
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Evolution of some orbital elements as a function of time for run 3b1 from Table 3: the eccentricity of the inner binary (*e*, left), and the mutual inclination between the two orbits (*i*', right).</figcaption></figure>

comment share a similar behavior, so that this constitutes a first classification of the runs.

### 2.2.1. Reference runs

What we hereafter refer to as *reference runs* are characterized by $q' \lesssim 800$ AU for Orbit BA, a moderate relative inclination ($i' \lesssim 40°$), and a current location of the GG Tau B binary more or less close to apoastron with respect to GG Tau A. These configurations are those which achieve the best fit of the shape of the out-coming circumbinary and the best secular stability of the stellar system. We describe a typical reference run, run 3b1 in Table 3. We first describe the orbital evolution of the stellar system over the time of the simulation. The system appears remarkably stable. Figure 2 shows the evolution of *e* and *i*' (recall that *i*' and $i_r$ coincide) as being the most critical parameters. These parameters exhibit a smooth oscillation of moderate amplitude around a mean value, and this is the most important variation we report here. The semi-major axes (not shown here) are extremely stable, but this will be the case for all the runs. It is well known that in non-resonant configurations (which is the case here), the semi-major axes are secular invariants. The outer orbit is very stable, with eccentricity change of less than 0.1%. We also have, as expected, a slow, regular precession of the periastra and lines of nodes of both orbits (not shown here), but this does not have much influence on the axisymmetric disk.

The evolution of the circumbinary disk in this environment is shown in Figs. 3 and 4. Figure 3 displays the number of remaining particles in the disk as a function of time. We see a rapid erosion of the disk in less than $2 \times 10^5$ yr, and a further stabilization of the remaining number of particles to ~10% of its initial number. This was expected and appears common to all runs. Choosing an initial disk significantly wider than the observed one, in particular with an outer edge sometimes above the periastron of the outer orbit, should inevitably lead to a considerable loss of particles. All lost particles appeared to have escaped from the system.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the number of remaining particles in the circumbinary disk as a function of time for run 3b1 from Table 3: ~10% of the particles initially present remain in the disk.</figcaption></figure>

Figure 4 describes the erosion of the initial disk, showing its radial profile at 4 different epochs. The remaining disk after erosion assumes a ring-like shape and closely resembles the observed one. We end up with a very sharp inner edge at ~180 AU (see also Fig. 6), i.e., exactly the same as what is observed, and a somewhat less sharp outer edge at ~300–350 AU. This again is a common feature with the observed disk, except that the real outer edge is rather ~260 AU than 350 AU. Interestingly, we see at $t = 1.5 \times 10^6$ yr that beyond the outer edge at 300 AU, there is still a tenuous tail of particles extending up to ~1000 AU. This tail disappeared by the end of the simulation. The $^{13}$CO $J = 2-1$ data (GDS99) imply the presence of such a tail and $1.5 \times 10^6$ yr could be a typical possible age for GG Tau. This suggests that the outer tail could be a transient structure due to the youth of the system.

This study first confirms the results of the two-body study (BD05). The inner edge of the disk appears where it is expected, given the inner orbit (Orbit AA2 here). This result is not affected by the presence of GG Tau B, nor by the resulting small amplitude libration of the eccentricity *e* of the inner orbit.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of the radial profile of the circumbinary disk (surface density of particles) as a function of time for run 3b1 from Table 3: the disk is rapidly eroded from the edges but a residual stable disk, close to the observed one, remains.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Semi-major axis – eccentricity profile of the circumbinary disk for run 3b1 from Table 3, at *t* = 1.5 Myr and *t* = 15 Myr. Each point corresponds to one particle remaining in the disk. At *t* = 1.5 Myr, a lot of very eccentric particles are still present beyond the outer edge of the disk, but they are ejected afterwards.</figcaption></figure>

The discrepancy between the inner edge of the disk and the orbital fit of GG Tau A remains unexplained. We discuss some possibilities for solving this problem. From a phenomenological point of view, the inner edge of the disk is obviously sculpted by the apoastron passages of the inner binary. Due to a small orbital period (360 yr), this sculpting is achieved quickly (in less than $10^4$ yr). The erosion of the outer edge is due to repeated periastron passages of the outer binary. Of course, as the orbital period is larger (55 775 yr), that sculpting takes a longer time and even if the outer edge at 300–350 AU appears after a few $10^5$ yr, its sculpting is not fully achieved before a few Myr. The present outer profile of the GG Tau disk may not yet have been achieved.

Finally, the ring-like disk we end up with corresponds to the largest one that can dynamically survive in such an environment. The orbital confinement of the disk can be seen in Fig. 5, which plots all remaining particles in the disk in a (semi-major axis – eccentricity) plane, at $t = 1.5 \times 10^6$ yr and at the end of the simulation. We see that in the semi-major axis, the outer edge of the disk is almost as sharp as the inner edge. The difference in the radial profile is due to the higher eccentricities

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Semi-major axis – inclination profile of the circumbinary disk for run 3b1 from Table 3, at *t* = 2 Myr.</figcaption></figure>

of the particles towards the outer edge than towards the inner edge. Thus, the outer edge is dynamically warmer than the inner edge. This fact shows up again in Fig. 6 which is similar to
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Fig. 4, but for runs 3b2 (**left**) and 3b3 (**right**) from Table 3, at *t* = 1.5 Myr (black) and *t* = 15 Myr (grey). The remaining disk is narrower than in Fig. 4 and almost disappears in run 3b2.</figcaption></figure>

Fig. 5, except that it plots now the disk in a (semi-major axis – inclination) plane, at *t* = 1.5 × 10⁷ yr. We clearly see the very sharp edges at 180 AU and 330 AU. Up to ~240 AU, the inclination distribution appears close to the initial one, but beyond that distance, the disk is somewhat excited by the outer binary.

The behavior reported here is common to all reference runs, i.e., there is always a residual disk in the expected semi-major axis range. The size of this disk depends however on the initial conditions of the run. Figure 7 shows the radial profile of the disk at *t* = 1.5 × 10⁶ yr and at the end of the simulation, but for runs 3b2 and 3b3 from Table 3. Compared to Fig. 4, the residual disk appears thinner. This is not surprising, as compared to run 3b1, runs 3b2 and 3b3 is characterized by a smaller periastron. In both cases, as for run 3b1, the orbital evolution of the system shows a remarkable stability; but due to a smaller periastron, the outer binary is deeper in the system, we logically expect the outer edge of the disk to be more eroded. This is confirmed by the simulation.

The residual disk obtained at *t* = 1.5 × 10⁶ yr in run 3b2 is very close to the observed one. The general shape is correct, and both edges almost exactly match the observations. This run is the one where we obtain the best fit of the observed disk. However, this disk is unstable and keeps being eroded a long time after. At the end of the simulation, it has almost disappeared, keeping only 100 particles from the 5000 initial ones (Fig. 7). Conversely, in run 3b3, characterized by a slightly larger periastron, the disk is somewhat thicker, but appears to better stabilize, suggesting a more realistic model.

For every run with “reference” initial conditions, we always end up with a residual disk comparable to the observed one, falling in the right semi-major axis range, with edges more or less close to the observations. The behavior reported here is thus very generic. Figure 8 shows the evolution of the radial profile of the disk in the same manner as Figs. 4 and 7, but for run 3b4 from Table 3. This run differs from run 3b1 only by the use of Orbit AA1 for GG Tau A instead of Orbit AA2 (see Table 1). This run is still classified as reference as the secular evolution of the whole system is qualitatively the same as for run 3b1, and because the outer profile of the disk is roughly identical to what is obtained with run 3b1 (compare to Fig. 4).

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 7, but for run 3b4 from Table 3, at *t* = 1.5 Myr and *t* = 15 Myr. The disk extends down to ~100 AU.</figcaption></figure>

However, as could be expected, the inner edge of the disk falls now at ~90 AU instead of 180 AU. This is obviously due to the smaller semi-major axis of Orbit AA1, and it is in perfect agreement with the study we performed in BD05 where only the GG Tau A binary was taken into account. This shows that the discrepancy we noted in BD05 between the observed inner edge of the disk and the best orbit deduced for the astrometric fit of GG Tau A (i.e., Orbit A1) is still present; this result is not affected by the perturbing action of the outer binary GG Tau B.

Similarly, if we change the outer orbit (BA), the resulting outer profile of the disk is affected, and may no longer match the observations. This occurs in Fig. 9, which is equivalent to Fig. 8, but for run 3b7 from Table 3. This run is a “large *q*’” run, characterized with respect to run 3b7 by a larger periastron (*q*') value for Orbit BA, and consequently a smaller eccentricity in order to fit the presently observed projected distance between the two binaries. Since in this case the two binaries have more distant approaches, the disk is less eroded from the outside. Its outer edge now falls at ~500–600 AU, which obviously does not match the observations. This behavior is common to all runs of the same kind, so that these orbital conditions must be rejected.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as Figs. 7 and 8, but for run 3b4 from Table 3, at *t* = 1.5 Myr and *t* = 15 Myr. The disk extends up to ~600 AU.</figcaption></figure>

### 2.2.2. The Kozai resonance

The runs described above indicate that the outer edge of the circumbinary disk is well simulated if we assume a periastron 600 AU $\lesssim q' \lesssim 800$ AU for the outer orbit. The inner edge location is compatible with Orbit AA2 for GG Tau A, but not with Orbit AA1. These conclusions hold if the orbital configuration of the two orbits does not undergo drastic secular changes. This was actually the case in the runs described above (see Fig. 2). All these runs assumed a moderate relative inclination between the two orbital planes. We now investigate the case of high initial *i*' values, as the resulting outcome may be drastically different from the previous case.

We start with run 3b5 (Tables 2 and 3), which mainly differs from runs 3b1–4 in the high initial *i'* value (80°). The main difference in the outcome concerns the orbital evolution of the system. Figure 10 shows the evolution of *e* (inner orbit) and of the mutual inclination *i*<sub>r</sub> as a function of time for run 3b5. It can be compared to Fig. 2 which is the same for run 3b1. Contrary to that case, we note here drastic changes to these parameters as a function of time. The eccentricity of the inner orbit now undergoes large amplitude oscillations over a period of ~8 × 10<sup>6</sup> yr, that bring it from ~0.35 to ~0.9. At the same time when *e* ≈ 0.9, the mutual inclination *i*<sub>r</sub> drops to ~42°.

This behavior is known as the *Kozai Resonance*. It is characteristic for triple systems with high inclinations, although it was described by Kozai (1962) for highly inclined comets. Under the effect of secular planetary perturbations (mainly from Jupiter), their orbit is subject to an evolution that drives it periodically to lower inclination but very high eccentricity, the semi-major axis remaining constant. The combined evolution of inclination and eccentricity is done in such a way that the action $\Theta = \sqrt{a(1-e^2)} \cos i$ is secularly preserved. This periodic evolution is characterized by librations of the argument of perihelion $\omega$ around 90° or -90°. As pointed out by Bailey et al. (1992), this mechanism is responsible for the origin of most sun-grazer comets in our Solar System, in particular those of the Kreutz group.

It is now well known that this dynamical mechanism is active in triple systems (see appendix). Harrington (1968) first

showed that it should be expected in hierarchical three-body system with high mutual inclination; Söderhjelm (1982) furthermore gave approximate analytical solutions for this problem, showing that the resonance should be active at $i_r \ge 40°$. It was recently reinvestigated in more details by Krymolowski & Mazeh (1999) and Ford et al. (2000). In Beust (2003), we used it as a tool to test the HJS integrator. The presence of the Kozai resonance is expected to affect the remaining circumbinary disk. Figure 11 shows the radial profile of the disk at three epochs, one before the eccentricity peak of the binary that occurs at $t = 2.7 \times 10^6$ yr, one after the peak, and at the end of the run. To a first approximation, the profiles appear close to the corresponding ones of Fig. 4 where no Kozai resonance was reported. The disk is not destroyed by the Kozai resonance. In the Kozai resonance, the outer orbit is only marginally affected by the dynamical evolution, so that we expect only minor consequences for the outer edge of the disk. Conversely, when the inner orbit undergoes a large eccentricity increase, we may expect the inner edge of the disk to be significantly eroded. Indeed, we note in Fig. 11 a significant erosion of the inner edge of the disk between $t = 1.5 \times 10^6$ yr and $t = 5 \times 10^6$ yr that is generally not reported in reference runs. This erosion is clearly related to the eccentricity peak of the inner orbit at $t = 2.7 \times 10^6$ yr. The effect is however not very large. The eccentricity of the inner orbit undergoes large amplitude variations due to the Kozai resonance, but its semi-major axis remains unchanged (this is a secular invariant in non-resonant problems). In BD05, we show that while the inner truncation of a circumbinary disk is a function of both parameters, it depends much more critically on the semi-major axis than on the eccentricity of the central binary.

This explains why the effect of the Kozai resonance on the profile of the circumbinary disk is rather limited. However, we note that after the eccentricity peak, the inner edge of the disk is pushed beyond 200 AU, which no longer matches the observations. This is nevertheless not enough to rule out the present run, because this additional truncation occurs after a time that is larger than the present age of the system.

A much more constraining output is the vertical shape of the disk. Due to the rapid precession of the orbital plane of GG Tau A, and to the direct perturbation by the highly inclined body GG Tau B, the disk achieves a much thicker vertical structure than in reference runs. This shows up in Table 3, where we note that the inclination dispersion $\bar{i}_d$ reaches 20° for run 3b5, which overcomes the 10° dictated by the observations. This is more visually illustrated in Fig. 12. In this figure, the vertical $(r, z)$ integrated shape of the disk is represented for runs 3b1 and 3b5, compared to the observational scale height law of Eq. (1). While the global shape for run 3b1 remains within the observational bounds, the disk for run 3b5 assumes a much more open and thicker structure, which does obviously not match the observation. Note that Fig. 12 is taken at $t = 2 \times 10^6$ yr for both runs, i.e., *before* the first eccentricity peak in run 3b5. For this reason, we think that the Kozai resonance scenario does probably not match the real situation.

This behavior is not specific to run 3b5; indeed, the Kozai resonance appears as soon as $i' \ge 40°$. Figures 13 and 14 describe the evolution of the system in the same manner as
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Same as Fig. 2, but for run 3b5 from Table 3: the eccentricity of the inner binary (e, left), and the mutual inclination between the two orbits (i<sub>r</sub>, right) as a function of time. Contrary to Fig. 2, e and i<sub>r</sub> exhibit drastic changes.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** The radial profile of the remaining circumbinary disk at *t* = 1.5 × 10<sup>6</sup> yr (black), *t* = 5 × 10<sup>6</sup> yr (dark grey) and *t* = 1.5 × 10<sup>7</sup> yr (light grey), for run 3b5 from Table 3. We note an inward erosion of the disk after a few Myr.</figcaption></figure>

Figs. 10 and 11 for run 3b5, but for run 3b6. Run 3b6 differs from run 3b5 only in the use of Orbit AA1 for GG Tau A instead of Orbit AA2 (Table 3). As for run 3b4, we expect for run 3b6 the inner edge of the disk to fall at ~100 AU rather than 180 AU, which is confirmed by Fig. 14. We also note in Fig. 13 a secular behavior typical of the Kozai resonance, except that the period of the Kozai modulation is now significantly longer than for run 3b5. We detect in Fig. 14 as in Fig. 11 an erosion of the inner edge of the disk due to the eccentricity peak.

The Kozai resonance could solve the paradox concerning the orbit of GG Tau A and the size of the inner gap. If the whole system was trapped in a Kozai resonance, then the present orbit of GG Tau A could be moderately eccentric (and correspond to the fit of the astrometric data, i.e., Orbit AA1), but the secular evolution could bring it periodically to high eccentricity values, leading to an inner erosion of the circumbinary disk that could not be achieved with a fixed orbit. Unfortunately, this does not work. There are three reasons for this. The first one appears in Figs. 14 and 11. Even if the inner edge of the disk is eroded

somewhat further, thanks to the Kozai resonance, the erosion remains moderate. This is due to the secular invariance of the semi-major axes.

Another reason is related to the vertical shape of the disk. As seen from Fig. 12, the resulting disk in the presence of the Kozai resonance is too thick to fulfill the observational constraints, and it can be seen from Table 3 that this is also true for run 3b6. This applies to all situations where the Kozai resonance is active.

Finally, we need to compare the age of the GG Tauri system with the period of the Kozai secular evolution. The age of the system does not exceed ~2 Myr. In Figs. 10 and 13, the period of the Kozai resonance is at least several Myr. This is far above the present age of the system. Thus, even if it is actually trapped in a Kozai resonance, with a present low *e* value, the GG Tau system is very unlikely to have undergone a high eccentricity peak in the last 2 Myr.

Runs 3b5 and 3b6 exhibit significantly different Kozai secular evolution periods. There is however no hope to find any other orbital configuration that would lead to a significantly smaller period (≤2 Myr). The period of the secular evolution cannot be expressed in a closed form, but theoretical orders of magnitude can be derived. A first order estimate of the period of the Kozai cycle may be expressed as (see appendix)

$$P_{\text{Kozai}} = \frac{M}{m_{\text{GG Tau B}}} \times \frac{P'^2}{P} \times (1 - e'^2)^{3/2} \times f, \quad (3)$$

where *M* is the total mass of the system, *m*<sub>GG Tau B</sub> is the mass of the outer body GG Tau B, *P* and *P'* are respectively the orbital periods of the inner and outer orbits, and *f* is a numerical factor of order unity that depends on the amplitude of the cycle.

This result is illustrated by the comparison between Figs. 10 and 13. When choosing Orbit AA1 for GG Tau A (run 3b6) instead of Orbit AA2 (run 3b5), we decrease *P* and increase *P*<sub>Kozai</sub> by a factor (62/32.4)<sup>3/2</sup> ≈ 2.65, which is roughly confirmed by Figs. 10 and 13. Once the inner orbit is fixed (e.g., Orbit AA1) the only way to reduce *P*<sub>Kozai</sub> is to choose the outer orbit in such a way to render the factor *P*'<sup>2</sup>(1 - *e*'<sup>2</sup>)<sup>3/2</sup> as small as possible, or similarly the term
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** The vertical $(r, z)$ profile of the remaining circumbinary disk at $t = 2 \times 10^6$ yr for runs 3b1 (**left**) and 3b5 (**right**). The thick grey lines show the observationally deduced scale height law (1) for comparison.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Same as Fig. 10, but for run 3b6 from Table 3. The behavior is similar but the period of the Kozai modulation is larger.</figcaption></figure>

$q'^3 ((1 + e')/(1 - e'))^{3/2}$. The outer orbit also needs to be compatible with the present projected separation and with the location of the outer edge of the disk. This fixes $q'$ to 700–800 AU. Then, the smallest $P_{\text{Kozai}}$ is obtained with the smallest possible $e'$, which is achieved saying that the outer orbit lies in the plane of the sky ($z' = 0$) and that it is currently at apoastron ($M' = 180^\circ$). Even doing this, various tests showed that we reduce $P_{\text{Kozai}}$ to 12–13 Myr, which is still far above the age of GG Tau.

### 2.2.3. Retrograde runs

Runs 3b8 and 3b9 from Table 3 describe “retrograde” runs, i.e., runs with the outer orbit almost retrograde with respect to that of GG Tau A. For what concerns our orbital parameters, this corresponds to $i' \gtrsim 140^\circ$.

We do not display the orbital evolution of the system for these runs, as it appears very similar to the prograde case (Fig. 2). Run 3b7 is identical to run 3b1, except that the inclination $i'$ was set to 160°, while run 3b9 corresponds similarly to run 3b2.

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Same as Fig. 11, but for run 3b6 from Table 3.</figcaption></figure>

The resulting disk profiles are shown in Figs. 15 and 16, which must be compared to Figs. 4 and 7. We see that in both cases, the disk is much less eroded outwards than in the corresponding prograde run. The residual ring of Fig. 15 (Run 3b8) is far too large to match the observations, and extends much
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** The radial profile of the circumbinary disk at *t* = 1.5 Myr and *t* = 15 Myr for run 3b8 from Table 3 (retrograde run). The disk is thicker than in corresponding prograde runs.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 17.** The vertical (*r*, *z*) profile of the disk like in Fig. 12, but for run 3b9 from Table 3. The disk is significantly thicker than observed.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Same as Fig. 15, but for run 3b9 from Table 3.</figcaption></figure>

further out than the disk in Fig. 4 (Run 3b1). Conversely, thanks to a smaller periastron value ($q' = 600$ AU), the residual disk in Fig. 16 (Run 3b9) appears closer to the observations. We nevertheless note that this disk is larger than the corresponding one of Fig. 7 (Run 3b2), showing that retrograde outer orbits actually less efficiently erode the disk than prograde orbits.

Another difference between prograde and retrograde runs appears in Table 3: retrograde runs are characterized by a larger inclination dispersion $i_d$ (20°), not compatible with the observational 10° limit. We illustrate this in Fig. 17 which shows the vertical profile of the disk as in Fig. 12, but for run 3b9. Obviously the disk is too thick to match the observations.

It thus seems that the outer edge of the disk is more stable but vertically thicker in the case of a retrograde orbit. Dynamically speaking, the problem of a test particle orbiting GG Tau A in the outer part of the disk can be viewed as close to a restricted, three-body problem, where GG Tau A is considered as a point mass. There is computational evidence that in the latter problem, situations where the particle orbits one of the two primaries in a retrograde manner are more stable than similar prograde configurations (Hénon 1970; Toomre & Toomre 1972; Wiegert & Holman 1997). The same kind of

stability was also found in hierarchical triple systems (Donnison & Mikulskis 1994).

The factor that tends to stabilize retrograde orbits and destabilize prograde ones in the restricted three-body problem is that the largest orbital disturbance takes place at conjunction with the perturbing body, i.e., when the particle passes between the two primaries. Consider now the motion of the particle in a rotating frame where the primaries have fixed positions. If the orbit is retrograde, then the relative velocity of the particle in that frame, with respect to the primaries, is higher. Hence the encounter time is much longer in prograde orbits than in retrograde ones. Moreover, when the orbit of the perturbing body (GG Tau B here) is eccentric, its instantaneous angular velocity with respect to the primary may match that of disk particles orbiting closer to the primary but on circular orbits. Thus we have a temporary, but highly destructive resonance, as noted by Toomre & Toomre (1972). This situation obviously does not concern retrograde particles which are consequently more stable.

Of course, this description holds only for situations close to instability, when the particle orbits close to its limiting Hill sphere. This concerns in our case the particles orbiting in the outer part of the circumbinary disk.

## 2.3. Non coplanar disks

As mentioned in BD05 and above, if we drop the assumption of coplanarity between the GG Tau A binary and the circumbinary disk, it is possible to find orbits for GG Tau A that are compatible with the astrometric data and with the inner edge location of the disk. The four solutions are listed in Table 1. We now explore the three-body dynamics with these solutions.

We tested many initial orbital configurations that we list in Table 4. For each run, we mention the orbital solution of Orbit AA (out of Table 1) that we use. Each run corresponds to a specific choice for Orbit BA, compatible with the present location of GG Tau B. For all these runs, we take $q' = 700$ AU and $M = 150°$. Then Orbit BA is fully specified giving its inclinations $i'$ and $i_r$ with respect to the plane of the circumbinary
<table><caption>**Table 4.** List of three-body runs with an initial non-coplanar disk described in this paper. For each run, we list the initial Orbit AA taken from Table 1, and the initial inclinations of Orbit AA $i'$ with respect to the disk and $i_r$ with respect to Orbit AA. Then we give as in Table 3 the approximate resulting boundaries of the main surviving ring, the approximate tilt angle between the midplane of the disk and the orbital plane of the central binary, the mean inclination $\bar{i}_d$ of the particles with respect to that midplane, and a key letter to indicate specific behaviors: “K” stands for Kozai (Kozai resonance), “R” for Retrograde (retrograde disk), “U” for Unstable disk (No disk left at the end of the run).</caption><thead><tr><th rowspan="2">Run<br>#</th><th rowspan="2">Orbit<br>AA</th><th rowspan="2">i'</th><th rowspan="2">i<sub>r</sub></th><th colspan="2">Ring boundaries (AU)</th><th rowspan="2">Tilt<br>angle</th><th rowspan="2">i&#773;<sub>d</sub></th><th rowspan="2">Remark</th></tr><tr><th>t = 2 Myr</th><th>t = 15 Myr</th></tr></thead><tbody><tr><td>3b10</td><td>AA4</td><td>0°</td><td>21.4°</td><td>180-320</td><td>190-280</td><td>≤10°</td><td>30°</td><td></td></tr><tr><td>3b11</td><td>AA4</td><td>20°</td><td>20°</td><td>180-330</td><td>190-280</td><td>13°</td><td>30°</td><td></td></tr><tr><td>3b12</td><td>AA4</td><td>80°</td><td>80°</td><td>180-350</td><td>200-300</td><td>~5°</td><td>40° K</td><td></td></tr><tr><td>3b13</td><td>AA4</td><td>160°</td><td>160°</td><td>180-750</td><td>180-600</td><td>~3°</td><td>32° R</td><td></td></tr><tr><td>3b14</td><td>AA4</td><td>20°</td><td>20°</td><td>180-330</td><td>190-290</td><td>12°</td><td>33°</td><td></td></tr><tr><td>3b15</td><td>AA4</td><td>80°</td><td>80°</td><td>180-450</td><td>170-400</td><td>3°</td><td>35°</td><td></td></tr><tr><td>3b16</td><td>AA4</td><td>160°</td><td>160°</td><td>180-650</td><td>180-600</td><td>~2°</td><td>31° R</td><td></td></tr><tr><td>3b17</td><td>AA3</td><td>0°</td><td>90.3°</td><td>170-350</td><td>180-250</td><td>~90°</td><td>80° K+U</td><td></td></tr><tr><td>3b18</td><td>AA3</td><td>20°</td><td>100°</td><td>170-400</td><td>180-300</td><td>~90°</td><td>80° K+U</td><td></td></tr><tr><td>3b19</td><td>AA3</td><td>20°</td><td>100°</td><td>170-400</td><td>180-290</td><td>~90°</td><td>80° K+U</td><td></td></tr><tr><td>3b20</td><td>AA3</td><td>80°</td><td>20°</td><td>160-280</td><td>170-240</td><td>~90°</td><td>80°</td><td></td></tr><tr><td>3b21</td><td>AA3</td><td>80°</td><td>20°</td><td>160-260</td><td>170-240</td><td>~90°</td><td>85°</td><td></td></tr><tr><td>3b22</td><td>AA3</td><td>80°</td><td>80°</td><td>170-500</td><td>170-340</td><td>~60°</td><td>85° K</td><td></td></tr><tr><td>3b23</td><td>AA3</td><td>80°</td><td>80°</td><td>160-320</td><td>170-260</td><td>~90°</td><td>85° K+U</td><td></td></tr><tr><td>3b24</td><td>AA3</td><td>80°</td><td>160°</td><td>170-380</td><td>170-270</td><td>~80°</td><td>90°</td><td></td></tr><tr><td>3b25</td><td>AA3</td><td>160°</td><td>80°</td><td>170-600</td><td>170-270</td><td>~100°</td><td>90° K+R+U</td><td></td></tr><tr><td>3b26</td><td>AA3</td><td>160°</td><td>80°</td><td>160-600</td><td>170-300</td><td>~80°</td><td>80° K+R+U</td><td></td></tr><tr><td>3b27</td><td>AA5</td><td>0°</td><td>23.7°</td><td>180-340</td><td>190-290</td><td>~12°</td><td>35°</td><td></td></tr><tr><td>3b28</td><td>AA5</td><td>20°</td><td>20°</td><td>180-360</td><td>180-280</td><td>~15°</td><td>40°</td><td></td></tr><tr><td>3b29</td><td>AA5</td><td>20°</td><td>20°</td><td>180-350</td><td>190-300</td><td>~15°</td><td>35°</td><td></td></tr><tr><td>3b30</td><td>AA5</td><td>160°</td><td>160°</td><td>180-650</td><td>180-600</td><td>~5°</td><td>35° R</td><td></td></tr><tr><td>3b31</td><td>AA5</td><td>160°</td><td>160°</td><td>180-700</td><td>180-550</td><td>~2°</td><td>35° R</td><td></td></tr><tr><td>3b32</td><td>AA5</td><td>80°</td><td>80°</td><td>180-450</td><td>200-340</td><td>~20°</td><td>80° K</td><td></td></tr><tr><td>3b33</td><td>AA5</td><td>80°</td><td>80°</td><td>190-400</td><td>200-300</td><td>~20°</td><td>70° K</td><td></td></tr></tbody></table>

disk and to Orbit AA respectively. For each run, we list the corresponding choice for $i'$ and $i_r$ in Table 4.

Of course not all combinations of $i'$ and $i_r$ are possible. Consider for instance run 3b10, which assumes Orbit AA4 and $i' = 0$, i.e., an outer orbit coplanar with the disk. Then the only possible choice for $i_r$ is the inclination between the disk and Orbit AA that is specific to Orbit AA4, i.e. $21.4°$. Similarly consider runs 3b17 to 3b26, which make use of Orbit AA3. That orbit is characterized by an inclination of $90.3°$ between the disk and Orbit AA. It is then not possible for example to find a solution for Orbit BA with $i' = 20°$ and $i_r = 20°$. The runs listed in Table 4 describe all typical possible configurations.

For any compatible choice of $i'$ and $i_r$, generally two orbital solutions for Orbit BA are found. In some cases the two solutions are very close, so that we test only one of them. But often the two solutions have very different semi-major axes and eccentricities, so that test both. This is why some lines in Table 4 share the same list of parameters.

For each run, we list in Table 4 (like in Table 3) the boundaries of the resulting rings of particles at $t = 2 \times 10^6$ yr and $t = 1.5 \times 10^7$ yr, the outcoming tilt angle between the ring and Orbit AA, and the inclination dispersion $\bar{i}_d$ of the particle disk with respect to its midplane. We finally give some key letters indicating the rough behavior: runs quoted with “K” are characterized by a pronounced Kozai resonance between the two orbits, leading in particular to a high amplitude modulation of the eccentricity $e$ of Orbit AA; “R” means a retrograde behavior, i.e., a less eroded disk thanks to a retrograde orbit of the GG Tau B perturber with respect to the disk; “U” stands for “Unstable disk”, i.e., a situation where virtually all particles in the disk are ejected from the system. The criterion is the following: a run is given the key letter “U” when less than 100 particles (among 50 000 initially) are left at the end of the simulation. Note that contrary to the coplanar case, here some runs may present more than one of these characteristics. Consider for instance runs 3b25 and 3b26. In these two cases, Orbit BA is initially almost coplanar with the disk, but retrograde ($i' = 160°$). We thus expect a retrograde-like behavior which is confirmed by the outer edge at 600 AU at $t = 2 \times 10^6$ yr. But as with Orbit AA3, the initial disk is almost perpendicular to the orbital plane of GG Tau A; the two orbits turn out to be necessarily almost perpendicular. Here we choose $i_r = 80°$. Consequently, we expect a strong Kozai resonance, which is confirmed by the simulation. The combination of these two characteristics leads to an instability of the disk after the first Kozai cycle. We note the same outcome if Orbit BA is initially prograde relative to the disk (Runs 3b17–3b19).

Finally, we find Kozai-like and Retrograde-like behaviors where we would expect them. It is more interesting to look at the evolution of the tilt angle between the ring and Orbit AA and the inclination dispersion $\bar{i}_d$. Recall that with Orbits AA4 and AA5 (Table 1) the disk is initially inclined by ~20° with respect to the orbital plane of GG Tau A, while in Orbits AA3 and AA6 it is nearly perpendicular to it. We note that for almost all runs involving Orbits AA4 and AA5, the resulting disk tends to realign with Orbit AA, but in the same time the inclination dispersion of the disk reaches fairly high values around ~30°. This may be explained as follows: all particles orbit initially the GG Tau A binary with inclinations around 20°. Under secular perturbations, the inclination is subject to fluctuations around this value, which may cause it to reach ~30° or more for some particles. Meanwhile, the line of nodes of all particle orbits is subject to precession, but the precession rate will depend on the orbit (basically its semi-major axis) of each particle. Therefore, if all particles share initially the same line of nodes (as they are part of an inclined disk), after some time their lines of nodes are randomly distributed, due to different precession rates. The net result is that the planar structure of the initial inclined disk is lost. The resulting disk assumes an open-cone shape coplanar with Orbit AA, but with a high opening angle (~30°). Finally this situation cannot match the observations, as the inclined disk appears not to keep its inclinations. Moreover, the 30° inclination dispersion is much higher than required.

If we consider now the runs with Orbits AA3 and AA6, we note that here the disk of particles remains perpendicular to the orbital plane of GG Tau A. In some cases, the disk is
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Views of the whole system (stars + particles) at $t = 4.5 \times 10^6$ yr for run 3b23 from Table 4. In the left plot, the orbit of GG Tau A is seen edge-on and the disk of particles is viewed face-on. In the right plot, the orbit of GG Tau A is viewed from the top and the disk is seen edge-on. The grey curve is the projection of Orbit BA.</figcaption></figure>

destroyed by tidal interactions (runs labelled “U”), but in other cases it survives. Hence a disk perpendicular to the orbital plane of GG Tau A could be stable. However, in all runs, the disk achieves a very high inclination dispersion (always ≥80°). In fact, with such a high $i_d$ value, can we still speak of a disk or just of a spherical shell of particles? In some cases we hardly recognize a disk, but in general, a ring-like structure is still visible. This is seen in Fig. 18, which shows two views of the spatial particle and body distribution for run 3b23 at $t = 4.5 \times 10^6$ yr. In one view the disk is seen face-on and in the other it is seen edge-on. We clearly identify a ring-like structure, although it can be seen in the edge-on view that is a very thick one. Thus, could this picture apply to the present status of the GG Tau system? Obviously no, because the disk, even if still detectable, is far too thick to be in agreement with the observational constraints.

Finally, all configurations involving disks that are not coplanar with the orbit of the GG Tau A binary must be rejected, because in all cases, whenever a disk survives, it turns out to be much too thick compared to the observed one. For the same reason, it seems that retrograde and Kozai situations need to be rejected. Therefore, it turns out that the most probable actual situation for the GG Tau A corresponds to our “reference” runs (runs 3b1–3b4 from Table 3), i.e., a disk coplanar with Orbit AA with a moderately inclined, prograde outer orbit. The periastron of Orbit BA must fall in the range 600–800 AU. The agreement between the astrometric measurements of Orbit AA and the location of the inner edge of the disk is still an open question.

# 3. The four body numerical study

## 3.1. Stability analysis

We now consider GG Tau B as a binary, by reintroducing the two bodies constituting it. We will assume the individual masses given by White et al. (1999), i.e., 0.12 and $0.044 M_\odot$ for GG Tau Ba and Bb respectively. As noted by

White et al. (1999), GG Tau Bb is a sub-stellar object. Considering GG Tau B as a binary means adding a third orbit to the system, i.e., the relative orbits between the two components of GG Tau B. For here on, we shall refer to this new orbit as Orbit BB, with doubly primed quantities ($a''$, $e''$...).

As for Orbit BA, we here again have only a projected distance ($1.5 \pm 0.08''$, i.e., 210 AU with $d = 140$ pc), and a position angle of $132 \pm 2.2^\circ$ (White & Ghez 2001). Four free parameters are needed to fully characterize that orbit. These will be its current mean anomaly $M''$, its eccentricity $e''$, an unknown altitude $z''$ of GG Tau Bb with respect to GG Tau Ba along the line of sight, and its inclination $i''$ with respect to the plane of the wide orbit BA between the two binaries.

We expect that the results of the three-body study, i.e., where GG Tau B was taken as a single body, should not be strongly affected by the four body study. Actually the most critical point concerning the four body study turns out to be the stability of the newly introduced orbit. GG Tau B is a much wider pair than GG Tau A. Its projected separation is ~200 AU. This means that its orbital semi-major axis $a''$ is at least 100–200 AU. The three-body study showed that the periastron $q'$ of Orbit BA should fall in the range 600–800 AU in order to allow a correct sculpting of the outer edge of the circumbinary disk. Thus the ratio $q'/a''$ is less than 10.

The ratio $q'/a''$ is expected to be the critical parameter controlling the stability of the system. For the stability of GG Tau B itself, the quadruple system may be dynamically viewed as a triple system constituting a single massive body (GG Tau A) orbiting a “central”, less massive binary (GG Tau B). There have been many studies about the stability of hierarchical triple systems. Harrington (1975) was the first to claim that the critical parameter was $q_{out}/a_{in}$, the ratio of the periastron of the outer orbit to the semi-major axis of the inner orbit. (i.e., exactly $q'/a''$ in the present case). He gave an empirical analytical stability criterion, $q_{out}/a_{in}$ required to be above a critical value to ensure stability. This question was reinvestigated by Donnison & Mikulskis (1994, 1995), and more recently by
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Summary of the stability study of the four-body GG Tau system: the lifetime of the GG Tau B pair is plotted as a function of the initial inclination *i''* assumed, for various orbital configurations. Each of the four frames corresponds to the use for Orbits AA and BA of one of the 4 “reference” 3-body configurations listed in Table 3. The different curves on each plot correspond to different combinations of initial eccentricity *e''* (0.2 or 0.7), and solution (1 or 2). The grey bar across the plots indicates the error bar for the age of the GG Tau system. Note that in some cases and for *i''* values close to 0° or 180°, no lifetime values are given because there is no orbital solution compatible with the present position of GG Tau B.</figcaption></figure>

Brasser (2002), on the basis of refined numerical studies, showing the limitations of the Harrington (1975) criterion, in particular in highly eccentric systems. All these studies fix the critical $q_{out}/a_{in}$ ratio between ~3 and ~10, depending on the relative masses of the individual components, and on the initial eccentricities and inclinations of the orbits. This shows that the GG Tau B pair should be close to instability.

We carried out a preliminary study where only the stability of the four-body system is investigated. We keep the four massive bodies but we consider no test particles. The circumbinary disk will be reintroduced afterwards. This allows to perform a large number of integrations very quickly with the HJS software. As we know that the system is close to instability, we reduce the time-step to 2 yr, i.e., 10 times less than in the three body study. In some of the following examples, we tested whenever instability was found that it was not a numerical artefact by i) again reducing the time step; and ii) performing a similar (but longer) integration using a conventional Burlish & Stoer (Press et al. 1992) integration of the same configuration. In all cases the instability was found at the same point of the simulation, showing that it is real (this appears not to be the

case if we keep the 20 yr time-step). We perform each integration as long as the GG Tau B pair remains bound ($e'' < 1$) and note the lifetime of the pair.

According to the $q'/a''$ criterion, the best stability should be achieved for any fixed $q'$ with the smallest $a''$ value possible. Small $a''$ values are obtained assuming that the orbit lies in the plane of the sky ($z'' = 0$), and that GG Tau B is currently close to apoastron ($M'' = 180°$). Moreover, high $e''$ values will subsequently lead to smaller $a''$ values. We thus decided to fix the initial values $M'' = 180°$ and $z'' = 0$ for all our tests. Indeed, for any integration, changing the initial $M''$ and $z''$ values should lead to less stability. For $e''$ we decided to systematically test for each possible combination a high ($e'' = 0.7$) and a low ($e'' = 0.2$) value.

The result of this study is summarized in Fig. 19. Each test is characterized by an initial choice for the three orbits. Each plot in Fig. 19 corresponds to a fixed use of a given solution for Orbits AA and BA, taken from the “reference” configurations listed in Table 3 (3b1–3b4). For each of these cases, various initial conditions are tested, depending on the choice of Orbit BB. The lifetime of that Orbit is plotted as a function of the initial
inclination $i''$. For each $i''$ value, we perform 4 different runs testing all the following combinations: $e''$ may be 0.2 or 0.7, and for each of these choices, there are still two possible orbital solutions, 1 and 2, for the GG Tau B orbit that are compatible with the present projected position. In each case, solution 1 is always the one leading to the larger inclination value with respect to the plane of the circumbinary disk.

The first conclusion is that the dynamical lifetime of the system is subject to considerable changes, by almost four orders of magnitude, depending on the initial conditions. In fact, some of the configurations tested may be indefinitely stable, but all integrations were stopped after 150 Myr. Conversely, some other configurations lead to an almost immediate instability. We consider that any configuration with a lifetime less than or comparable to the age of GG Tau is unrealistic.

The plots in Fig. 19 are very noisy. We performed integrations starting with all possible integer values for $i''$ in degrees between 0 and 180°. Between two adjacent runs with initial $i''$ values separated by only 1°, the outcomes lifetime may be very different. This reveals the highly chaotic nature of the evolution of Orbit BB. However, despite this intrinsic chaos, some general trends appear. First, the configurations with initial $e'' = 0.7$ (thick lines) appear in general more stable than those with $e'' = 0.2$. This is in agreement with our preliminary analysis of the $q'/a''$ criterion.

The same criterion also explains why some orbital configurations of Orbits AA and BA lead to better stability, irrespective of the choice of Orbit BB. Obviously, the choice of the 3b2 configuration almost never leads to stability (Fig. 19). Indeed, this configuration assumes the smallest periastron ($q' = 600$ AU) of our sample. The 3b3 configuration ($q' = 700$ AU) leads to better stability and the 3b1 and 3b4 ones ($q' = 800$ AU) to even better stability.

Finally, the stability of the system appears to be highly dependent on the initial inclination $i''$ choice. Clearly, retrograde orbits ($i'' \ge 100°$) are generally more stable than prograde orbits corresponding to the same configuration. This relies to the same phenomenon as we noted in the previous section: in the same way as the orbits of the particles in the circumbinary disk were more stable when they were retrograde with respect to Orbit BA, here Orbit BB is more stable when it is retrograde with respect to that Orbit too. We have some more or less isolated configurations that are stable compared to the age of the GG Tau system despite a prograde Orbit BB ($i'' < 90°$). But most of the stable configurations fall in the retrograde regime.

We note that the results with the 3b1 and 3b4 configurations are very close. These configurations differ only the solution assumed for Orbit AA (AA1 or AA2). This illustrates the fact that, as seen from GG Tau B, GG Tau A mainly acts as a single perturbing body.

## 3.2. Full study

From this study we see that the dynamical stability of the outer pair GG Tau B is in itself a very constraining parameter for the whole system. For instance, we may from now reject the 3b2 configuration characterized by $q' = 600$ AU, because there

<table><caption>**Table 5.** The initial configurations taken for Orbit BB in the full four-body runs described in Table 6. The inclination $i''$ is given relative to the wide orbit BA.</caption><thead><tr><th>Orbit designation</th><th>Semi-major axis $a''$</th><th>Eccentricity $e''$</th><th>Inclination $i''$</th><th>Altitude $z''$</th><th>Current mean anomaly $M''$</th></tr></thead><tbody><tr><td>BB1</td><td>122 AU</td><td>0.7</td><td>151°</td><td>0</td><td>180°</td></tr><tr><td>BB2</td><td>173 AU</td><td>0.2</td><td>168°</td><td>0</td><td>180°</td></tr><tr><td>BB3</td><td>173 AU</td><td>0.2</td><td>149°</td><td>0</td><td>180°</td></tr></tbody></table>

<table><caption>**Table 6.** List of four-body runs performed with a particle disk. The parameters listed for each run are identical to those listed for the three-body runs listed in Table 3. Each run corresponds to the use of one of the reference configurations listed in Table 3 to which a retrograde stable Orbit BB is added. The comparison three-body run is given in the last column of the table.</caption><thead><tr><th>Run #</th><th>Orbit combination</th><th>Ring boundaries (AU)<br>t = 2 Myr<br>t = 15 Myr</th><th>Tilt angle</th><th>$i_d$</th><th>Comp. 3-body run</th></tr></thead><tbody><tr><td>4b1</td><td>AA2+BA1+BB1</td><td>180–400</td><td>180–330</td><td>≤1°</td><td>8° 3b1</td></tr><tr><td>4b2</td><td>AA2+BA3+BB2</td><td>180–340</td><td>190–280</td><td>≤5°</td><td>17° 3b3</td></tr><tr><td>4b3</td><td>AA1+BA1+BB3</td><td>90–400</td><td>90–350</td><td>≤1°</td><td>12° 3b4</td></tr></tbody></table>

is no corresponding stable solution for Orbit BB. We must have now $q' \ge 700$ AU to ensure stability for Orbit BB. Moreover, from Fig. 19 Orbit BB is very likely to be retrograde with respect to Orbit BA. Therefore, the most probable orbital configuration for the whole GG Tau system should be the following: Reference conditions for Orbit BA, i.e, a moderate $i'$ inclination and 700 AU $\le q' \le 800$ AU; an Orbit BB retrograde with respect to Orbit BA, with a high $e''$ eccentricity, and currently at apoastron. We now present a full study (with disk) of this configuration.

We performed three typical four-body simulations, each one corresponding to the use of one of the reference three-body configurations listed in Table 3, to which we add a retrograde choice for Orbit BB that exhibits stability in the preliminary study over much more than the 15 Myr time-span of the full simulations. These orbits are listed in Table 5 and the result of the runs is summarized in Table 6, which may be directly compared to Table 3.

Details about these runs are shown in Figs. 20–23. Figure 20 shows the secular evolution of Orbit BB. The evolution concerning run 4b3 is not shown as it is identical to that of run 4b2. The large eccentricity variations reveal the chaotic nature of the dynamics. Even the semi-major axis exhibits some chaotic behavior. However it remains close to its initial value, showing that the orbit is nevertheless stable, though chaotic.

Figures 21 and 22 show the evolution of the radial profile of the disk for the three runs of Table 6, in the same manner as Figs. 4, 7 and 8 which show the same for the corresponding three-body runs. In all cases, the radial profile of the disk is close to that obtained in the corresponding three-body run. This is confirmed by the data of Table 6, which shows in comparison to Fig. 3 that the inclination and inclination dispersion parameters of the disk are very similar. More specifically,
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Orbital evolution (semi-major axis and eccentricity) of Orbit BB in runs 4b1 and 4b2 from Table 6. The evolution is significantly chaotic though stable.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 21.** The radial profile of the circumbinary disk at *t* = 1.5 Myr (black) and *t* = 15 Myr (grey), for run 4b1.</figcaption></figure>

<caption>**Fig. 22.** Same as Fig. 21, but for run 4b2.</caption>

the outer edge of the disk is in general slightly more eroded
in the four-body runs than in the corresponding three-body
run. This is particularly true for run 4b2 (Fig. 22) compared
to run 3b3 (Fig. 23). At *t* = 15 Myr, the disk in run 4b2 is much
less crowded (though having a similar radial extension) than in
run 3b3. This effect is much less pronounced for the other runs,
and this reveals the origin of the phenomenon. The additional

erosion of the disk is due to the successive inwards excursions
of the components of GG Tau B towards the outer edge of the
disk, their center of mass still following the given Orbit BA
common to the corresponding three-body and four-body runs.
The effect is more important for run 4b2 because this run as-
sumes a periastron *q*' = 700 AU for Orbit BA compared to
*q*' = 800 AU for the other runs.
<figure><img src="image_23.png" /><figcaption>Fig. 23. Same as Fig. 21, but for run 4b3.</figcaption></figure>

In our runs, we derive disks with profiles that are fairly close to the observed one. This concerns the disk at *t* = 1.5 Myr for run 4b1 (Fig. 21) and the disk at *t* = 15 Myr for run 4b2 (Fig. 22). The inclination parameters of these disks (Table 6) are also compatible with the observations: the disks are coplanar with Orbit AA and the inclination dispersion *i*<sub>d</sub> remains acceptable (though somewhat larger than the required 10° for run 3b2). The orbital configuration quoted here, common to the three runs described, is able to sculpt a circumbinary disk with characteristics that match the observations. Of course, the disk obtained with run 4b3 is not compatible with the data, but this is due to the use of Orbit AA1 (Table 1) and concerns only the inner edge of the disk. This was already true in the three-body runs and illustrates the discrepancy between the observed inner edge of the disk and the astrometric measurements of the GG Tau A pair. Comparing runs 4b1 and 4b3 shows that changing the orbit of GG Tau A (Orbit AA) has virtually no influence on the outer profile of the disk. Only the inner edge is affected. There is therefore no hope to derive additional constraints on Orbit AA via the general dynamical behavior of the whole system. As noted above, the same conclusion is derived from the stability analysis of Orbit BB (Fig. 19).

There is also another interesting byproduct of this work: why is there a circumbinary disk around GG Tau A and not around GG Tau B? In order to investigate this issue, we tried to add a second disk of particles orbiting GG Tau B, starting with the same conditions as run 4b1. We do not display the corresponding results as they are identical to those of run 4b1. The additional disk of particles appears to be fully destroyed within less than 10⁵ yr. A circumbinary disk orbiting GG Tau B is thus not stable. This is not surprising, as GG Tau B itself is only marginally stable. Note that the fact that no circumbinary disk is to be expected around GG Tau B does not exclude that circumstellar material could be present orbiting the individual components of GG Tau B. White et al. (1999) note spectroscopic indications for accretion onto both components of GG Tau B. This is another issue that could be investigated in future work.

4. Discussion

The most important outcome of the preceding numerical study is that taking into account the constraints arising from the circumbinary disk and from the stability analysis of GG Tau B, it is possible to derive a likely configuration for the whole system that is able to sculpt a disk matching the observations:

1. The circumbinary disk is coplanar with GG Tau A (Orbit AA).

2. The periastron *q*′ of the wide Orbit BA must fall in the range 700–800 AU.

3. The Orbit BA is more or less coplanar with GG Tau A and the disk, the relative inclination *i*<sub>r</sub> being in any case <40° in order to avoid the Kozai resonance.

4. The orbit of the GG Tau B pair (Orbit BB) is fairly chaotic and very probably retrograde with respect to the two other orbits; it is also likely not to be significantly inclined with respect to the plane of the sky, in order to ensure better stability.

Nearly all simulations presented in this paper show that the sculpting process of the disk is not completed at the suspected age of the system (1–2 Myr), as the disk appears significantly further eroded at the end of the simulations at *t* = 15 Myr. This suggests that the circumbinary disk presently observed around GG Tau A could be a transient feature that could undergo a significant further erosion in the few next Myr. However, viscous processes in the disk, not taken into account in our pure *N*-body simulations, could help stabilize a residual disk against further erosion, as they would tend to reduce the velocity dispersion within the disk. Obviously a further study, possibly SPH (Smoothed Particle Hydrodynamics) is needed to take into account viscosity. The global shape of the stable disks we derive with our *N*-body calculations should not be drastically affected by gas drag.

As shown in Beust (2003), symplectic integration (done with HJS) and SPH simulation performed by Artymowicz & Lubow (1994) of the sculpting process of standard circumbinary disks yield similar results. This shows that the sculpting is basically gravitationally driven; gas dissipation acts only as a minor correction. Even if the dynamical picture we describe here is more complex, this basic result should hold. The motivation for using a symplectic integrator for GG Tauri in the present study is that the whole system is characterized by several very different dynamical time-scales. The sculpting process of the inner edge of the disk is much more rapid than that of the outer part. Hence if we want to resolve both we need to adopt a time-step dictated by the smaller time-scale, and integrate over a time span dictated by the longer. This is easily achieved with symplectic integration but still hardly reachable with SPH.

The dynamical stability of the whole four-body system is not ensured over a long time-span. The chaotic nature of Orbit BB shows that the GG Tau B pair could be tidally disrupted within a few Myr. At least one of the components of GG Tau B would then escape the system. A good way to constrain this problem, and to confirm the retrograde motion of Orbit BB, would be to gather some astrometric measurements
of GG Tau B. This is however much more difficult than for GG Tau A. Given the above quoted orbital characteristics (Table 5) the period of Orbit BB should fall in the range 3000–6000 yr. Assuming that this orbit lies close to the plane of the sky, the rate of change of the relative position angle of the components of GG Tau B should not exceed ~0.1°/yr. It could easily be one order of magnitude less if the orbit is eccentric and currently at apoastron. Given the error bar on the measurement of that position angle (132 ± 2.2° White & Ghez 2001), we would have to wait several decades before detecting any change. This is nevertheless not unreachable, especially if the accuracy of the astrometric measurement is improved in the near future. Therefore initiating an astrometric monitoring of the GG Tau B pair, as already done for several years for GG Tau A, could be appropriate. We could obtain valuable information about the whole system, including the disk orbiting GG Tau A, by gathering measurements on both pairs, not only GG Tau A.

The major discrepancy between the orbit we fit for GG Tau A from the astrometric data (Orbit AA1) and the one needed to account for the location of the inner edge of the disk (Orbit AA2) remains. Further monitoring of the pair is needed to specify this fact and in particular to check the accuracy of error bars on the astrometric data. However, if the present astrometric solution (Orbit AA1) is confirmed in the near future, then the discrepancy will remain unsolved. In BD05, we discuss some ways to solve the discrepancy. We note first that invoking viscosity in the disk does not solve the problem, as it renders it even more drastic. We suggest that the past secular evolution of Orbit AA could possibly explain the discrepancy, e.g. by an inward migration of the orbit. However, such a migration would also imply a similar migration of the inner edge of the disk, so that the problem remains unsolved.

Moreover from the present study we see that Orbit AA is remarkably stable (Fig. 2). Even with the eccentricity modulation induced by the Kozai resonance, there is no way to force the erosion of the disk up to 180 AU instead of 90 AU starting from Orbit AA1.

Finally, in BD05, we also suggest that the disk could be non-coplanar with GG Tau A. The present dynamical study shows that this possibility must be ruled out, as the resulting disks do not match the observations.

## 5. Conclusion

Thanks to the abundance of available data and the flexibility of our HJS code, we are able to constrain the dynamics of the whole GG Tauri A and B system which appears to be hierarchical. This is the first attempt of such a detailed dynamical study. Even if some points remain unclear, our study allows us to derive a more likely scenario that is able to explain the inner and outer ring shape of the circumbinary disk. Our analysis of the dynamical stability also reveals that the CB ring is almost coplanar with the GG Tau A system.

Surprisingly, the study of the outer pair GG Tau B is of interest due to their marginal stability, as it provides constraints

on the whole system. Its orbit is likely retrograde. Moreover, the presence of a circumbinary disk around GG Tau B is very unlikely and gravitational instabilities would destroy it very quickly.

Finally, the discrepancy between the inner edge of the disk and the existing astrometric data of GG Tau A remains and should be investigated when more astrometric data become available. More astrometric data will also allow us to refine the determination of the orbit of GG Tau B.

*Acknowledgements.* We thank Arnaud Pierens for communication of results prior to publication. All the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI).

## References

Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651

Bailey, M. E., Chambers, J. E., & Hahn, G. 1992, A&A, 257, 315

Beust, H. 2003, A&A, 400, 1129

Beust, H., & Dutrey, A. 2005, A&A, 439, 585 (BD05)

Beust, H., Corporon, P., Siess, L., Forestini, M., & Lagrange, A.-M. 1997, A&A, 320, 478

Brasser, R. 2002, MNRAS, 332, 723

Chambers, J. E. 1999, MNRAS, 304, 793

Donnison, J. R., & Mikulskis, D. F. 1994, MNRAS, 266, 25

Donnison, J. R., & Mikulskis, D. F. 1995, MNRAS, 272, 1

Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067

Dutrey, A., Guilloteau, S., & Simon, M. 1994, A&A, 286, 149

Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385

Guilloteau, S., Dutrey, A., & Simon, M. 1999, A&A, 348, 570 (GDS99)

Harrington, R. S. 1968, AJ, 73, 190

Harrington, R. S. 1975, AJ, 80, 1081

Hénon, M. 1970, A&A, 9, 24

Kawabe, R., Ishiguro, M., Omodaka, T., Kitamura, Y., & Miyama, S. M. 1993, ApJ, 404, L63

Kinoshita, H., & Nakai, H. 1999, Celest. Mech., 75, 125

Kozai, Y. 1962, AJ, 67, 591

Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720

Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18

McCabe, C., Duchêne, G., & Ghez, A. M. 2002, ApJ, 575, 974

Pierens, A., Huré, J.-M., Dutrey, A., Guilloteau, S. 2005, in preparation

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes (Cambridge Univ. Press.)

Roddier, C., Roddier, F., Northcott, M. J., Graves, J. E., & Jim, K. 1996, ApJ, 463, 326

Söderhjelm, S. 1982, A&A, 107, 54

Tamazia, V. S., Docobo, J. A., White, R. J., & Woitas, J. 2002, ApJ, 578, 925

Toomre, A., & Toomre, J. 1972, ApJ, 178, 623

White, R. J., & Ghez, A. M. 2001, ApJ, 556, 265

White, R. J., Ghez, A. M., Reid, N. I., & Schultz, G. 1999, ApJ, 520, 811

Wiegert, P. A., & Holman, M. J. 1997, AJ, 113, 1445

Wisdom, J., & Holman, M. 1991, AJ, 102, 1528
**Online Material**
## Appendix A: The Kozai resonance in the stellar three body system

Our stellar three-body system is fully described by the description of the two orbits (AA and BA). Let us call **r** and **r'** their radius vectors. Starting from its expression in barycentric coordinates, a straightforward algebra gives the expression of the Hamiltonian *H* of the system (Beust et al. 1997):

$$H = -\frac{Gm_0m_1}{2a} - \frac{G(m_0 + m_1)m_2}{2a'} + Gm_0m_2\left(\frac{1}{r'} - \frac{1}{r_{02}}\right) + Gm_1m_2\left(\frac{1}{r'} - \frac{1}{r_{12}}\right), \quad (A.1)$$

where *a* and *a'* are the semi-major axes of the orbits, *G* is the gravitational constant, *m*₀ and *m*₁ are the masses of the central bodies (GG Tau A) and *m*₂ is the mass of the outer body (GG Tau B); *r*₀₂ and *r*₁₂ are the distances from bodies 0 to 2 and 1 to 2 respectively. This is a peculiar application of the canonical transformation given in Beust (2003) that helps to express the Hamiltonian of the general *N*-body hierarchical problem into generalized Jacobi coordinates.

In the expression of *H*, *r*₀₁ and *r*₀₂ are implicitly functions of **r** and **r'**. More specifically, we get

$$r_{02} = r' \sqrt{1 + \frac{2m_1 \cos \beta}{m_0 + m_1} \frac{r}{r'} + \frac{m_1^2}{(m_0 + m_1)^2} \frac{r^2}{r'^2}}, \quad (A.2)$$

$$r_{12} = r' \sqrt{1 - \frac{2m_0 \cos \beta}{m_0 + m_1} \frac{r}{r'} + \frac{m_0^2}{(m_0 + m_1)^2} \frac{r^2}{r'^2}}, \quad (A.3)$$

where *β* is the angle between **r** and **r'**. The first two terms of *H* represent the uncorrelated Keplerian orbits while the other terms constitute the disturbing function. As a function of *r*, *r'* and cos*β*, *H* contains short period terms. The Kozai Hamiltonian itself is the double time average of the full Hamiltonian *H*. It describes accurately the secular dynamics of the three-body system as long as the two orbits are not locked in resonance, which is the case here:

$$\bar{H} = \frac{1}{4\pi^2} \int_0^{2\pi} \int_0^{2\pi} H(l, l') \, dl \, dl', \quad (A.4)$$

where *l* and *l'* are the mean anomalies. In fact, deriving conceptually $\bar{H}$ is achieved through a canonical transformation called the von Zeipel transformation, which leads to an Hamiltonian totally free of short period terms. The orbital elements appearing in $\bar{H}$ may be seen as time-averaged versions of the oscillating ones. From a technical point of view, $\bar{H}$ is just obtained from Eq. (A.4).

Unfortunately, $\bar{H}$ cannot be expressed analytically. The classical procedure is to expand the Hamiltonian in ascending powers of *α* = *a*/*a*', taking advantage from the fact that *α* ≪ 1, or equivalently *r*/*r*' ≪ 1. From the above expressions of *r*₀₂ and *r*₁₂, we see that 1/*r*₀₂ and 1/*r*₁₂ naturally expand in ascending powers of *r*'/*r* using Legendre Polynomials. The final result for *H* is:

$$H = -\frac{Gm_0m_1}{2a} - \frac{G(m_0 + m_1)m_2}{2a'} - \frac{Gm_0m_1m_2}{a'} \times \left\{ \sum_{j=2}^{+\infty} \alpha^j \frac{m_0^{j-1}(-m_1)^{j-1}}{(m_0+m_1)^j} \left(\frac{r}{a}\right)^j \left(\frac{a'}{r'}\right)^{j+1} P_j(\cos\beta) \right\}, \quad (A.5)$$

where the *P*<sub>*j*</sub>'s is the *j*th Legendre polynomial. This expression of the Hamiltonian is given in Harrington (1968), Krymolowski & Mazeh (1999) and Ford et al. (2000).

Using *α* ≪ 1, the expansion is then truncated up to some given order and averaged over the mean longitudes. Theoretically, once it is truncated to some given finite order, the expansion of *H* can be averaged analytically. However, the complexity of the expression of the truncated $\bar{H}$ increases dramatically with the order of the expansion. The expansion limited to second order only (*j* = 2) is usually called quadrupolar while the expansion up to *j* = 3 is called octupolar.

In order to perform the integration efficiently, it is better to describe the orbits relative to the natural referential frame of the system, where the OZ axis is parallel to the constant, total angular momentum of the system. It is also convenient to introduce the mutual inclination *i*<sub>r</sub> between the two orbits. In this referential frame, it is well known that the longitude of nodes Ω and Ω' of the two orbits verify Ω - Ω' = π, which yields more simple expressions. The resulting quadrupolar expansion is

$$\bar{H}_{\text{quad}} = -\frac{1}{16} \frac{Gm_0 m_2 m_1}{m_0 + m_1} \frac{\alpha^2}{a' \sqrt{1 - e'^2}} \times [(2 + 3e^2)(3 \cos^2 i_r - 1) + 15e^2 \sin^2 i_r \cos 2\omega]. \quad (A.6)$$

This Hamiltonian is given by Ford et al. (2000) and Krymolowski & Mazeh (1999). Kinoshita & Nakai (1999) also derive a similar approximation of the Kozai Hamiltonian for the asteroid case. Ford et al. (2000) and Krymolowski & Mazeh (1999) also give full expressions for the octupolar expansion. Interestingly, Krymolowski & Mazeh (1999) retain additional terms ∝ *α*⁷/² in the octupolar expansion that arise from the von Zeipel transformation, but according to Ford et al. (2000), these terms have a negligible role in the dynamics. No such terms are present at the quadrupolar expansion level.

Interestingly, $\bar{H}_{quad}$ does not depend on *ω*', which shows that limited to quadrupolar expansion, the eccentricity *e*' of the outer orbit is a secular invariant. The variations of *ω* and *i*<sub>r</sub> are coupled with those of *e* by conservation of the total angular momentum of the system. Hence *i*<sub>r</sub> can be eliminated from the expression of $\bar{H}_{quad}$, so that $\bar{H}_{quad}$ turns out to have only one degree of freedom (it is thus integrable), coupling the secular variations of *e* and *ω*. The dynamics of the system can be deduced from level curves of $\bar{H}_{quad}$ in a (*a*, *e*) plane for a given value of the angular momentum. The Kozai resonance itself is characterized by a libration regime of *ω* around ±π/2, coupled with high amplitude oscillations of *e*.

We wish to only derive an estimate of the period of the Kozai cycle, so that we will limit ourselves to the quadrupolar expansion. In order to derive the equations of motion, the Hamiltonian $\bar{H}_{quad}$ should be given as a function of the conjugate Delaunay elements of both orbits. This is done in Ford et al. (2000) and Krymolowski & Mazeh (1999). The variation rate of *e* reads

$$\begin{aligned} \frac{de}{dt} &= -\frac{\sqrt{1 - e^2}}{e \sqrt{a(m_0 + m_1)}} \frac{\partial \bar{H}_{\text{quad}}}{\partial \omega} \\ &= -\frac{15}{8} \frac{m_2 \sqrt{G}}{\sqrt{m_0 + m_1}} \frac{a^{3/2} e \sqrt{1 - e^2}}{a^{\prime 3} (1 - e^{\prime 2})^{3/2}} \sin^2 i_r \sin 2\omega \end{aligned} \quad (A.7)$$
$$ = -\frac{15\pi}{4} \frac{m_2}{m_0+m_1+m_2} \frac{e\sqrt{1-e^2}}{(1-e'^2)^{3/2}} \frac{P}{P'^2} \sin^2 i_r \sin 2\omega, $$

where P and P' are respectively the orbital periods of the inner and outer orbits. The period of the Kozai modulation reads

$$ P_{\text{Kozai}} = \frac{8}{15\pi} \frac{m_0 + m_1 + m_2}{m_2} \frac{P'^2}{P} (1 - e'^2)^{3/2} \times \left| \int_{e_{\text{min}}}^{e_{\text{max}}} \frac{\mathrm{d}e}{\sin^2 i_r \sin 2\omega e \sqrt{1 - e^2}} \right|. \quad (A.8) $$

$e_{\text{min}}$ and $e_{\text{max}}$ are the minimum and maximum eccentricity values reached in the Kozai cycle. The integral appearing in that expression cannot be expressed in closed form; it can nevertheless be obtained using the Weierstrass's zeta function (Kozai 1962). This integral is related to the amplitude of the Kozai

eccentricity modulation and is most cases of order unity. For a given amplitude (e.g., $e_{\text{min}} = 0.3$ and $e_{\text{max}} = 0.95$), we should not expect it to vary drastically if we change the initial orbital configuration. Finally, $P_{\text{Kozai}}$ turns out to depend mainly on the factor $(P'^2/P)(1 - e'^2)^{3/2}$, confirming the result of Söderhjelm (1982) who gave $P'^2/P$ as a typical time-scale for that dynamics.

In the quadrupolar approximation, the eccentricity oscillations are strictly periodic, but in the full system, the higher order terms induce both amplitude and period fluctuations in the successive cycles (Ford et al. 2000; Krymolowski & Mazeh 1999); the order of magnitude $P_{\text{Kozai}}$ remains nevertheless the same. All these high order effects are very accurately reproduced with the HJS integrator, as shown in Beust (2003).