A&A 387, 550-559 (2002)
DOI: 10.1051/0004-6361:20020407
R. Günther - W. Kley
Institut für Astronomie & Astrophysik, Abt. Computational Physics, Auf der Morgenstelle 10, 72076 Tübingen, Germany
Received 21 January 2002 / Accepted 14 March 2002
Abstract
We study the evolution of circumbinary disks surrounding
classical T Tau stars.
High resolution numerical simulations are employed to model
a system consisting of a central eccentric binary star within
an accretion disk. The disk is assumed to be infinitesimally thin,
however a detailed energy balance including viscous heating
and radiative cooling is applied.
A novel numerical approach using a parallelized Dual-Grid technique
on two different coordinate systems has been implemented.
Physical parameters of the setup are chosen to model
the close system of DQ Tau and AK Sco,
as well as the wider systems
of GG Tau and UY Aur. Our main findings are
for the tight binaries a substantial flow of material through the
disk gap which is accreted onto the central stars in
a phase dependent process. We are able to constrain the parameters
of the systems by matching both accretion rates and derived spectral
energy distributions to observational data where available.
Key words: accretion, accretion disks - binaries: general - hydrodynamics - methods: numerical
Apart from many spectroscopic binary systems with surrounding disks like DQ Tau and AK Sco, there have been only two observations of directly imaged wide optical binary systems with circumbinary disks, namely GG Tau (Dutrey et al. 1994) and UY Aur (Duvert et al. 1998). The observations of the spectroscopic binaries usually include emission data from the infrared and optical part of the spectrum. Observations of the luminosities of the boundary layer are typically used to estimate accretion rates onto the binary stars (Gullbring et al. 1998; Hartigan et al. 1995).
Bate & Bonnell (1997) model the early phase of a molecular cloud collapsing onto a protobinary system with an SPH code and derive properties of the created circumbinary and circumstellar disks. We are interested in the following phase of the evolution when the envelope has been cleared. By modeling the evolution we compare fully developed circumbinary and circumstellar disks and their properties with observational data.
In a system consisting of a stellar binary surrounded by an accretion disk there is a continuous exchange of angular momentum between the binary and the disk. As the binary is orbiting with a shorter period than the disk material it is generating positive angular momentum waves in the disk which, upon dissipating, deposit their energy and angular momentum in the disk. Hence, there is a continuous transfer of angular momentum from the binary to the outer disk. Upon receiving this angular momentum, the inner disk speeds up, and the material is receding from the binary, leaving a region of very low density around the central binary, a so called inner gap. At the same time, the loss of angular momentum leads to a secular shrinkage of the semi-major axis of the binary. The width and form of the gap depends on disk parameter such as temperature and viscosity and on the orbital parameter (eccentricity) of the binary as well (Artymowicz et al. 1991; Artymowicz & Lubow 1994).
First calculations dealing with the evolution of circumbinary disks so far have been using the numerical method of smoothed particle hydrodynamics. Applying such a scheme Artymowicz & Lubow (1996) could demonstrate that for sufficiently high viscosities and temperatures material from the disk is able to penetrate into the inner gap region of the disk to become eventually accreted by one of the binary stars. This process which is facilitated by a gas flow along the saddle points of the potential leads to a pulsed accretion whose magnitude and phase dependence for eccentric binaries was estimated. However, when using only a small number of particles all having the same mass only a limited resolution is achievable.
Later Rozyczka & Laughlin (1997) used a finite difference method in attacking this problem. Using a constant kinematic viscosity coefficient and eccentric binaries they also found, in agreement with Artymowicz & Lubow (1996), a pulsed accretion flow across the gap, with most of the accretion occurring onto the lower mass secondary star.
In this paper we extend these computations and solve explicitly a time dependent energy equation which includes the effects of viscous heating and radiative cooling. We aim at the structure and dynamics of the disk as well as the gas flow in the close vicinity of the binary star. To that purpose we utilize a newly developed method which enables us to cover the whole spatial domain. For the first time we perform long-time integration of the complete system covering several hundred orbital periods of the binary and compare the properties of the evolved systems with observational data such as spectral energy distributions in the infrared and optical bands and accretion rates estimated from luminosities.
In the next section we layout the physical model followed by some remarks concerning important numerical issues (Sect. 3). We describe the setup of the various numerical models in Sect. 4. The main results are presented in Sect. 5 and our conclusions are given in Sect. 6.
The gas in the disk is non-self-gravitating and is orbiting a binary system with stellar masses M_{1} and M_{2} which are modeled by the gravitational potential of two solid spheres. The binary is fixed on a Keplerian orbit with given semimajor axis a and eccentricity e. Hence, we neglect self-gravity and the gravitational backreaction of the disk onto the stars. This is justified because the total mass of the disk , within the simulated region which extends up to several tens of semi-major axis a of the binary, is typically much smaller than that of the central binary. Hence, the timescale for a change in a and e is indeed much longer than the timescales considered here (Artymowicz et al. 1991). For the majority of the computations we employ an inertial coordinate system, however test calculations were performed in a system corotating with the binary as well.
The gravitational potential ,
generated by the binary stars,
is given by
For the kinematic shear viscosity
we assume in this work
an -model (Shakura & Sunyaev 1973) in the form
(5) |
(6) |
(7) |
(10) |
(11) |
(12) |
(13) |
As the calculations have shown, the system behaves extremely dynamically with very strong, moving gradients at the inner edge of the disk. To model such a situation, the viscosity has to be treated implicitly to ensure numerical stability. The coupling of and requires a solution method similar to that in Kley (1989). For stability purposes an additional implicit artificial viscosity has to be included in the computations (Kley 1999).
This code has been employed successfully in related simulations dealing with embedded planets in disks (Kley 1999; Kley 2000). The numerical details of the present simulations are presented in Günther (2001).
In the main part the domain is covered by a cylindrical grid. On this (outer) grid the equations are formulated in polar coordinates, and the solution technique guarantees exact conservation of global angular momentum. The conservation property is absolutely crucial in obtaining physically reliable results for long term disk evolution (Kley 1998).
However, such a grid cannot be extended to the very center and
leaves a hole in the middle. For typical computations
this does not pose serious problems, but in this case we
are dealing with eccentric binary systems in the middle. Hence, the
stars or their circumstellar disks may easily cross the inner grid boundary causing strong
irregularities. In any case a central hole does not allow for an accurate
determination of the individual mass accretion rates onto or a possible
overflow between the stars.
Therefore, using a Cartesian grid centered on the origin enables an accurate
calculation of the mass flows in the vicinity of the stars.
Applying this technique, a high local
resolution can be obtained in the center as well.
An example of such a grid structure is displayed in Fig. 1.
Figure 1: Grid system of a typical calculation. Shown is the grid structure of the Cartesian and cylindrical grid at the center of the computational domain. Circumstellar disks of the system can be seen at the top and bottom edges of the picture. | |
Open with DEXTER |
Related numerical schemes using a nested grid technique, though not in different coordinate systems, have been applied in astrophysical simulations by a number of authors (eg. Ruffert 1992; Yorke et al. 1993).
The necessary equations (in Cartesian and cylindrical coordinates) are then integrated, independently and in parallel, on the two grids, and after each timestep the necessary information in the overlap region has to be exchanged. Restrictions are imposed on the time step for stability reasons, the Courant-Friedrichs-Lewy (CFL) condition must be fulfilled during each integration, on each grid.
Test calculations using a single Cartesian grid covering the whole domain have shown that the non-conservation of angular momentum yields serious deviations from known analytical results in particular for longer integration times.
As the grids possess a different geometry, exact conservation of momenta and energy cannot be guaranteed in the conversion/interpolation process which uses linear interpolation to exchange density, temperature and velocity. Exchanging density, internal energy and momentum was also implemented and leads to more correct conservation of momenta and energy. However, the advantage of a complete coverage of the physical space in between the two stars clearly outweighs this disadvantage.
The shock experiment initial setup is a step function of density
and temperature with the floor starting at r_{0} (about 10 grid cells outside of the inner boundary of the cylindrical grid) and covering the whole Cartesian grid:
Figure 2: Radial shock in first (infalling) phase right after crossing the grid boundaries. Surface density is color coded, the grid structure is shown. | |
Open with DEXTER |
Figure 3: Radial shock in second (explosion) phase right after crossing the grid boundaries again. Surface density is color coded, the imaged domain is the same as in Fig. 2. | |
Open with DEXTER |
From now on, we refer to non-dimensional units. All the lengths are expressed in units of the semi-major axis a of the binary. This is constant because the time scale on which the binary orbital elements change due to mass and angular momentum transfer is much longer than the dynamical time scales considered here (see Artymowicz & Lubow 2001). Masses are in units of the solar mass and time is given in units of the binary's orbital period.
The whole azimuthal range of the disk is taken into account by considering a computational domain represented by , where typically in units of the binary separation. The cylindrical grid covers a region with represented by up to meshpoints, and the central Cartesian grid has up to meshpoints covering with an additional small overlap.
The region R out of which matter is removed is given as a fraction (0.2) of the Roche lobe size for the wide systems and as the stellar cores for the tight spectroscopic systems. The accretion process typically involves only a few grid cells on each side of the stars for the wide systems and the whole resolved stellar core (up to grid-cells) for the tight systems, making it a locally confined process.
Two methods for deciding whether a grid cell is targeted for accreting
a part of its mass were implemented. Accreting a constant
fraction f per time of the mass exceeding a defined minimum density
leads to a very smooth accretion process.
The accreted mass per iteration is computed by
The inferred accretion rates in any event compare favorably with the average disk accretion rate, as estimated from the radial mass fluxes through different surfaces, and the evolution of the total mass found in the circumstellar disks.
Spectral energy distributions are generated decoupled from the
dynamic evolution of the disk as a post-processing step. Here the
model from Adams et al. (1988) is extended to work on
data generated by our numerical simulation. The flux at the observer
is obtained by integrating over the disk surface assuming local
black-body radiation
The integration is limited to regions outside of the stellar cores as temperatures and densities from inside the cores can be nothing but wrong. The star emission is accounted for by adding a black-body spectrum from an appropriate effective star temperature.
Figure 4: Initial gap profile composed out of and the gap function Eq. (18). The figure shows a gap at with . | |
Open with DEXTER |
The gravitational potential in Eq. (4) is phased-in smoothly from an initial central potential for the total mass, M_{1} + M_{2}, during an initialization phase lasting typically ten orbital periods. This allows the system to adjust smoothly from the semi-stationary setup of a single central star to the binary star situation.
For the boundary conditions, periodicity is imposed at and . We set outflowing boundary conditions at the outer radial border ( ), i.e. the radial velocity gradient at is zero, and truncated to positive values. Angular velocity there equals the unperturbed initial Keplerian value.
For single grid calculations with a cylindrical grid we use zero gradient boundary conditions for both velocity components and the density at the inner border.
In Fig. 5 you can see a circumbinary disk after 90 orbits evolution time in almost perfect quasi-stationary state. The deviation from line-symmetry is small. Only with very high resolution simulations we observe spiral features in the outer regions of the disk, while spiral arms originating from the inner gap of the circumbinary disk flowing onto the circumstellar disks can be observed even for low resolutions. Those spiral arms feeding the circumstellar disks remain stable for the whole period of a non-eccentric binary.
The accretion rates for both accretion mechanisms onto the two stars are nearly identical and do not show any significant variations for different numerical parameters. Thus, our implementation of the accretion mechanism is robust.
Figure 5: Evolution of an equal mass binary after 90 orbital periods. Color encoding is logarithmic surface density. Both the inner part of the circumbinary disk and the circumstellar disks are visible. | |
Open with DEXTER |
Changing model parameters such as ,
gap width, or disk size we
obtain characteristic variations of the resulting spectral energy
distribution of an evolved disk. Figs. 6-8
show dependencies of these parameters on the spectral energy distribution
that are used to fit parameters of the spectroscopic systems to match
observed spectra.
Figure 6: Dependency of the gap width on the SED. | |
Open with DEXTER |
Figure 7: Dependency of the disk size on the SED. | |
Open with DEXTER |
Figure 8: Dependency of on the SED. | |
Open with DEXTER |
AK Scorpii | DQ Tauri | |
P | ||
e | ||
q | ||
i | ||
a | ||
T_{*} | 6500 K | 4000 K |
R_{*} | ||
M_{*} | ||
d | ||
Figure 9: AK Sco circumbinary disk after 82.0 orbital periods in apastron. Color coding is , the size of the stars reflects the actual stellar radii, the length scales are in AU. | |
Open with DEXTER |
Hereafter we discuss models for the two spectroscopic
young binary stars DQ Tau and AK Sco which
have a separation of only a fraction of an .
The physical parameter for the systems are taken from
Andersen et al. (1989),
Favata et al. (1998), Jensen & Mathieu (1997)
for AK Sco and Mathieu et al. (1997)
for DQ Tau.
The relevant information for the
simulations has been summarized in Table 1,
disk masses
are integrated densities over the computational
domain specified by the disk radii .
Images showing the disk configuration in apastron and periastron
phase are shown in Figs. 9 and 10.
Figure 10: DQ Tau circumbinary disk after 85.5 orbital periods in periastron. Color coding is , the size of the stars reflects the actual stellar radii, the length scales are in AU. |
Hartigan | Gullbring | simulation | |
DQ Tau | |||
AK Sco | n.a. | n.a. | |
GG Tau | |||
UY Aur |
Figure 11: Time-dependent accretion rate of the AK Sco system. The orbital phase curve schematically shows the distance of the two stars. | |
Open with DEXTER |
Generally, for eccentric systems we observe spiral features in the outer circumbinary disk only for very high resolution simulations, while again spiral arms feed the circumstellar disks even with low resolutions. Opposed to the non-eccentric systems the spiral arms develop during apastron phase and are torn off together with the outer parts of the circumstellar disks after periastron phase for high eccentric systems. This tearing off is caused by raising centrifugal forces acting on the circumstellar material.
For both systems we compared the accretion rates resulting from the
simulation with accretion rates derived from of observations as done
by Hartigan et al. (1995) and Gullbring et al. (1998)
(see Table 2). As no data was available for AK Sco,
a comparison with observational data was not possible, but as both
spectroscopic systems have similar system parameters, accretion rates
of the same order of magnitude can be expected.
Figure 12: Time-dependent accretion rate of the DQ Tau system. The orbital phase curve schematically shows the distance of the two stars. | |
Open with DEXTER |
As can be seen from Figs. 11 and 12 the accretion process happens synchronized with periastron phase of the binaries. This is in agreement with numerical simulations from Artymowicz & Lubow (1996) and Rozyczka & Laughlin (1997) and with observations which show excess luminosity at these phases (Mathieu et al. 1997).
From our simulations we can derive an averaged accretion rate of for the primary, for the secondary of AK Sco and for the primary, for the secondary of DQ Tau. These accretion rates are in good agreement with the observational rates.
It is expected that the luminosity increases in periastron because tidal compression of the disk leads to a temperature rise with subsequent luminosity increase. We tried to generate luminosity curves to show that indeed excess luminosity is produced at periastron phase, but unfortunately those curves show both brightening and darkening dependent on the implemented accretion mechanism. At the very low averaged accretion rates of the disk is very optically thin. In that regime the implemented radiative loss mechanism is not very accurate and thus these results are not meaningful.
For the spectroscopic binaries we can generate spectral energy distributions out of the simulation data which can be made to match the observed ones by manually adjusting initial model parameters such as disk mass, gap width and density profile d. For DQ Tau Fig. 13 shows the spectral energy distribution composed out of the disk and the star, Fig. 14 shows the one for the AK Sco system.
The effective stellar temperatures used in the model have been adjusted to match the observational data. They differ from the quoted values in the literature in Table 1. This may be due to difficulties in separating individual stellar and disk contributions to the total flux in models and observations.
Also occultation is not accounted for on adding the star spectrum to the disk.
Figure 13: Spectral energy distribution of DQ Tau (d=2.0), crosses show observational data taken from Mathieu et al. (1997). | |
Open with DEXTER |
Figure 14: Spectral energy distribution of AK Sco (d=1.5), crosses show observational data taken from Jensen & Mathieu (1997). | |
Open with DEXTER |
UY Aurigae | GG Tauri | |
P | 2074 yr | 490 d |
e | 0.13 | 0.25 |
q | 0.63 | 0.77 |
i | ||
a | ||
T_{*} | 4000 K | 3800 K |
R_{*} | ||
M_{*} | ||
d | ||
Figure 15: UY Aur circumbinary disk after 20.5 orbital periods. Color coding is , the length scales are in AU. | |
Open with DEXTER |
From our simulations we derive averaged accretion rates for both UY Aur and GG Tau that are comparable with the observed quantities in Table 2. In the case of UY Aur we get for the primary and for the secondary, respectively. Simulations of GG Tau lead to accretion rates of for the primary and for the secondary. As expected from the difference of the disk masses ( vs. ) the accretion rate for GG Tau is one order of magnitude less than the one for UY Aur.
Accretion is reinforced at periastron phases as for the spectroscopic systems, but events are an order of magnitude lower for GG Tau and nearly vanish for UY Aur due to the lower eccentricity of the systems. Also both systems show more phase dependent accretion for the secondary than for the primary.
In Figs. 15 and 16 you can see
the inner parts of the simulated circumbinary disks. The circumstellar
disks remain intact over the whole period.
It was not possible to generate meaningful spectral
energy distributions of the disks as they are too cold for these wide
systems. Also observational data is lacking for these systems.
Figure 16: GG Tau circumbinary disk after 56.4 orbital periods. Color coding is , the length scales are in AU. | |
Open with DEXTER |
To obtain a detailed comparison with observed light curves in particular the phase dependence a more elaborate model of the radiative loss of the very thin material in the circumstellar disks is required. In that case 3d simulations may be necessary which are beyond present day computational resources.
Acknowledgements
Part of this work was funded by the German Science Foundation (DFG) under SFB 382 Simulation physikalischer Prozesse auf Höchstleistungsrechnern. We would like to thank Dr. Hubert Klahr for many stimulating discussions during this project.