Radiative transfer in circumstellar disks

Radiative transfer in circumstellar disks

I. 1D models for GQ Lupi
S. D. Hügelmeyer Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    S. Dreizler Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    P. H. Hauschildt Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany    A. Seifahrt Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    D. Homeier Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    T. Barman Lowell Observatory, 1400 W Mars Hill Rd, Flagstaff, AZ 86001, USA
Received 17 December 2008 / Accepted 10 March 2009
Key Words.:
Accretion, accretion disks – Radiative transfer
Key Words.:
Radiative transfer - Accretion, accretion disks - Methods: numerical - Techniques: spectroscopic

We present a new code for the calculation of the 1D structure and synthetic spectra of accretion disks. The code is an extension of the general purpose stellar atmosphere code PHOENIX and is therefore capable of including extensive lists of atomic and molecular lines as well as dust in the calculations. We assume that the average viscosity can be represented by a critical Reynolds number in a geometrically thin disk and solve the structure and radiative transfer equations for a number of disk rings in the vertical direction. The combination of these rings provides the total disk structure and spectrum. Since the warm inner regions of protoplanetary disks show a rich molecular spectrum, they are well suited for a spectral analysis with our models. In this paper we test our code by comparing our models with high-resolution VLT CRIRES spectra of the T Tauri star GQ Lup.

1 Introduction

Gas and dust disks are common objects; they can be observed around a variety of objects such as very young stars (e. g. T Tauri and Herbig Ae/Be stars), evolved binaries (cataclysmic variables), and even black holes. With the discovery of the first extrasolar planets over years ago, the interest in protoplanetary disks has increased. Disk properties such as density, temperature, and chemical composition effect the process of planet formation and therefore also the characteristics of the planets. For very young stars (ages of ), the accretion of disk matter onto the star is an internal source of energy for the disk. The mass of such a disk is typically a few percent of the stellar mass and mass accretion rates are of the order of . Furthermore, irradiation by the central star heats the outer layers of the disk which, in turn, radiate the impinging energy away. Therefore, the direct detection of protoplanetary disks provides the possibility to observe the earliest evolutionary stages of planet formation.

Tremendous efforts have been made to model the structure and more importantly for the comparison with observations radiative transfer in accretion disks. Kriz & Hubeny (1986), Shaviv & Wehrse (1986), Hubeny (1990) and others went beyond the vertically averaged models of Shakura & Syunyaev (1973) and Lynden-Bell & Pringle (1974) or models using the diffusion approximation (e. g.  Meyer & Meyer-Hofmeister 1982; Cannizzo et al. 1982) to obtain full numerical solutions for the structure of and the radiative transfer in accretion disks. Since then, these models have reached a high degree of sophistication (for an overview of protoplanetary disk models see Dullemond et al. 2007).

Even though instruments like the VLTI constitute a major improvement for the observation of spatially extended objects, most of our information about the physics of protoplanetary disks comes from spatially unresolved spectra. NIRSPEC observations of protoplanetary disks have shown the capabilities of IR spectroscopy (Najita et al. 2003). With the ESO IR spectrograph CRIRES (Kaeufl et al. 2004) such observations have become available to a larger community of astronomers. Dust continuum radiative transfer (RT) calculations (e. g. Wolf 2001) can reproduce the overall structure of the dust disk out to a few hundred AU. The more abundant gas in the disk, however, cannot be predicted by these models. The warm inner disks (temperatures between to a few ) provide a special laboratory to study the gas structure because temperatures and densities are adequate to produce molecular spectral lines visible in the near- to mid-infrared. High-resolution observations in combination with model spectra enable us to obtain kinematic information about the gas since line profiles are governed by the velocity field in the disk. Another interesting observation is the agreement of mean inner gas disk radii and orbital radii of short-period extrasolar planets (Carr 2007). To further investigate this and other phenomena, detailed gas and dust models of the warm inner disk regions are necessary.

Our paper is structured as follows: In Section 2 we will explain the construction of our model structure and synthetic spectra. We will concentrate on the basic concepts and highlight the implementation of the solution. In Section 3 we will present synthetic spectra for the disk of the T Tauri star GQ Lup and compare these with CRIRES infrared observations to show the potential of our 1D disk model code. Section 4 will give a summary and an outlook of our work.

2 Models

We have developed a circumstellar disk radiative transfer code as an extension of the well-tested model atmosphere package PHOENIX (Hauschildt & Baron 1999). PHOENIX can handle very large atomic and molecular line lists and blanketing due to several hundred million individual lines is treated in the direct opacity sampling method. Dust is included in the models presented here by treating condensate formation under the assumption of chemical equilibrium and phase-equilibrium for several hundred species. I. e. all dust monomers exceeding the local saturation pressure as defined by thermal equilibrium are allowed to condense (this implies a supersaturation ratio ; Allard et al. 2001, cf. also Helling et al. 2008 for a comparison of different condensation treatments).

The grain opacity is calculated from the most important refractory condensates using optical data for a total of 50 different species. Absorption and scattering are calculated in the Mie formalism, assuming a given particle size distribution for a mixture of pure spherical grains, following the general setup of the Dusty set of PHOENIX atmosphere models (Allard et al. 2001, a plot with relative abundances of the most important species in our disk models is shown in Fig. 9). This equilibrium assumption is a good approximation in the inner optically thick layers of the well-mixed disk atmosphere. In the low-density outer layers, non-thermal effects are more important; including these effects in the DISK version of PHOENIX is planned for the future. In the later phases of disk evolution grain growth will become important and may cause departures from a homogeneous size distribution. Partial pressures for the different atoms, molecules, and dust species are calculated with the new “Astrophysical Chemical Equilibrium Solver” (ACES) equation of state (Barman, in preparation) which allows us to reach temperature regimes as low as a few tens of degrees Kelvin.

In our models, the disk region considered is divided into rings (see Fig. 1) and a plane-parallel disk atmosphere is computed between the midplane and the top of the disk for a given number of quadrature points (Gaussian angles) for each ring independently. Here denotes the angle between the characteristic and the normal to the disk plane.

We adopt the standard accretion model for geometrically thin disks, i. e. the disk height is much smaller than the disk ring radius . This assumption decouples the treatment of vertical and radial disk structure because the vertical structure is in quasi-static equilibrium compared to time scales for the radial motion of gas. Matter is assumed to rotate with Keplerian velocities and viscous shear decelerates inner and accelerates outer parts leading to accretion of matter and outward transportation of angular momentum. Molecular viscosity is too small to provide the observed mass accretion rates. Thermal convection in accretion disks was investigated by different authors using various methods which are summarised in Klahr (2007). The results show that thermal convection is unlikely the dominant source of turbulence and even leads to inward transport of angular momentum. Furthermore, a heating source is required to drive the convection. Bell & Lin (1994) argue that convective instabilities can only occur at temperatures  K or , i. e. temperature regimes which our models do not reach (see Sect. 3). The magnetorotational instability (MRI) introduced by poloidal magnetic fields in weakly ionised disk matter (Balbus & Hawley 1991) is the currently favoured origin of viscosity but its effect on the thermal structure of the disk cannot be easily described or parametrised. Even though temperatures  K are necessary to thermally ionise disk material, cosmic ray ionisation is possible at surface densities (Gammie 1996) which is true for all our models presented below. The mean viscous dissipation is

Figure 1: Disk ring structure as adopted for our calculations. The radius of the rings increases exponentially. The left panel shows a face-on view of a disk, the right one a vertical cut viewed edge-on (height is not to scale). The dotted lines are the radii for which the models are calculated while the solid lines show the borders of a disk ring. The disk structure is assumed to be constant over the ring width.

often modelled as an “alpha-viscosity” resulting in a vertically-averaged viscosity


(Shakura & Syunyaev 1973) where is the angular momentum transfer efficiency, the sound speed, and the pressure scale height. Alternatively, a turbulent viscosity can be modelled assuming a critical Reynolds number


(Lynden-Bell & Pringle 1974). This second model has the advantage, that the calculation of the mean viscous dissipation is decoupled from the thermal structure of the disk and is adopted here. Both of these formalisms allow one to account for the effect of viscosity on the disk structure without the need to describe its origin in detail.

2.1 Start models

In order to start a new DISK disk calculation either an existing model can be used as input or a grey LTE start model is constructed. In the latter case we follow the approach of Hubeny (1990).

The model is initialised by setting up a mass depth scale , where is the number of layers, which will be kept fixed also for later calculations with the same input parameters. The scale is equally spaced on a logarithmic grid between the column mass at the midplane of the disk


and the outer value , which is an input parameter. This spacing is done for points, while for the remaining layers each point , , is positioned halfway between and , with . This is done to provide numerical stability in the iterative process because of the steep gradient near the midplane. A typical value for is or for .

In a next step, we calculate the depth-dependent viscosity


assuming a value of determined by an assumed constant critical Reynolds number (Eq. 2), where is a free parameter (Kriz & Hubeny 1986). The smaller the value of , the lower the temperature. However, if irradiation is noticeably strong, the effect of on the structure has negligible influence on the spectrum. Lynden-Bell & Pringle (1974) argue that the Reynolds number has to be chosen to equal the critical one for the onset of turbulence. We usually take which leads to .

We further assume an isothermal density structure – the sound speed associated with the gas pressure and the flux mean opacity are independent of height – which lets one derive a simple analytic expression for the density at each depth point . The relation


is then inverted to convert our mass depth variable into a height above the midplane of the disk .

Finally we employ the formal LTE solution derived in Hubeny (1990) and the Rosseland opacity tables of Ferguson et al. (2005) to get a temperature structure by iterating


with until the --structure is consistent for each layer. In Eq. 6 we set equal to and is the optical thickness at the midplane of the disk (see Sec. III b of Hubeny 1990).

2.2 Iterative procedure

After a restart model is read in or a start structure has been computed, the hydrostatic equation, the radiative transfer (RT) equation, and the energy balance equation have to be solved. This is done iteratively until convergence in flux is obtained. Our termination criterion is


where is the expected “mechanical” flux released by the disk and the current flux value. This convergence criterion in Eq. 7 is not always reached and for a sufficiently small the calculation is then stopped after a given number of iterations. Typical errors in are of the order of .

Figure 2: Temperature correction scheme for a not irradiated disk atmosphere with , , , , and . The top plot of the left panel shows the relative temperature change for each iteration in the top most layer. The bottom plot is the temperature in that layer. The right plot shows relative temperature corrections as a function of layer number for a chosen set of iteration steps (see legend). The dotted line shows the correction in the first and the dashed line in the last iteration. The correction is limited to in the first four iterations to avoid overcorrections. In this example convergence in the Unsöld-Lucy temperature correction is reached after iterations.

2.2.1 Hydrostatic equation

The hydrostatic equation is one of the basic equations that govern the structure of a stable disk. Since the mass of the disk is much smaller than that of the central object, we can neglect self-gravitation of the disk. Assuming that the radial component of the central star’s gravitation and the centrifugal forces of the rotating disk just cancel each other out, we obtain the vertical hydrostatic equation for a thin disk


where is the sum of gas pressure and radiation pressure . It is now convenient to use the relation and replace Eq. 8 by


This way, we eliminate the height (which is not known a-priori) and introduce the sound speed , which is taken from the previous iteration and only depends on the temperature . Therefore, the error in and cancel out and we obtain a more stable version of Eq. 8. This second-order differential equation is then decomposed into two coupled first-order differential equations. With the inner boundary condition (symmetry at the midplane) and the outer boundary condition derived in Hubeny (1990), we solve the two point boundary value problem with a simple shooting method.

2.2.2 Radiative transfer

We solve the plane-parallel radiative transfer equation


for a given number of Gaussian angles , , with corresponding quadrature weights from the midplane (, ) to the top of the disk (, ). The upper boundary condition for Eq. 10 is


and the lower boundary condition due to symmetry conditions is


In Eq. 11, the expression denotes the impinging radiation from the central star at the top of the disk atmosphere. The radiation field is computed using the Accelerated Lambda Iteration (ALI; see, e. g., Hauschildt & Baron 1999). Hence, we have to incorporate the boundary conditions given in Eqs. 11 and 12 into the formal solution (FS) of the transfer equation (10) and into the Approximate Lambda Operator (ALO) , which is introduced to improve the convergence rate of the Lambda Iteration by reducing the eigenvalues of the amplification matrix. The FS can be written as


where is the source function and is the thermal coupling parameter. The -operator is split according to


and we can rewrite Eq. 13 as


The non-local -operator is calculated according to Hauschildt et al. (1994).

2.2.3 Energy balance

The standard accretion disk model demands that all the energy that is produced by viscous dissipation, i. e. mechanical energy, is released in form of radiative and convective energy, viz.


For a viscous disk, the left hand term can be written as


(Kriz & Hubeny 1986) and the radiative energy is expressed by


We shall neglect the convective energy term from here on because radiation is dominating the temperature structure and convection seems to have little effect (D’Alessio et al. 1998). We then derive an Unsöld-Lucy class temperature correction scheme (Lucy 1964) by integrating the second frequency-integrated moment of the radiative transfer equation




The introduction of the Eddington factors


and insertion of the first frequency-integrated moment of the radiative transfer equation gives the temperature correction law




the absorption mean and Planck mean opacity. An example iteration history is shown in Fig. 2.

2.2.4 Irradiation

Irradiation by the central star plays an important role in the determination of the temperature and height profile of a protoplanetary accretion disk. Therefore, we have taken special care to treat this effect in detail. The impinging intensity onto the surface of the disk (see Eq. 11) is determined by first calculating the slope of the disk surface at ring radius by computing the height of the disk according to our start model (see Section 2.1) at with . For each Gaussian angle and corresponding integration weight the projected surface fraction on the star is determined. This area is then subdivided and a mean intensity reaching the disk surface is evaluated for a star with a blackbody spectrum according to


where is the number of surface segments for Gauss angle , is the distance from the area element to the disk surface, and the angle under which the radiation leaves the star. The factor is needed since the the irradiation onto the disk is unidirectional and not isotropic. Limb darkening is considered by the factor . Alternatively, a PHOENIX spectrum can be used as irradiation source. The irradiation geometry is shown in Fig. 3.

2.2.5 Spectrum

Having determined a set of self-consistent disk ring structures, we

Figure 3: Irradiation geometry as adopted for our calculations. We consider a star with radius at a distance from a disk ring with a height and a slope of the disk surface . For each Gaussian angle and its corresponding integration weight the projected surface fraction on the star is determined and the irradiation intensity is calculated. The angle under which radiation leaves the surface of the star is .

calculate a last iteration with a larger number of Gaussian angles (usually ) and a fine wavelength spacing (up to in the wavelength range of interest) to obtain a high-resolution model spectrum. Since the Keplerian rotation velocity and the inclination angle of the disk are not considered in the single ring spectra yet, we compute a combined spectrum according to


In Eq. 26 we assumed that the disk is axis-symmetric and that the intensity is constant for all radii between inner and outer

Figure 4: Optical depth structure for a disk ring model with , , , , and . The left plot shows the run of with height ( is the midplane of the disk). The asterisks denote the optical depth at the center of the line, the crosses are for a continuum point (see right panel). The dotted line at marks the region where the line and the continuum become optically thick, i. e. where the radiation that we see comes from. The plot on the right is a temperature--diagram showing at which temperatures and column masses the optical depth for the line (dashed line) and the continuum (dotted line) become unity.

radii for each ring. In addition to the integration over all disk rings, the influence of the disk’s rotation on the line profile is taken into account by applying the Doppler shift


to the line. Here the velocity is determined for a given inclination and for a set of disk ring segments with azimuthal angle . We use disk ring segments, i. e.  steps in , to determine the rotationally broadened line profile. This method is a simple way to determine line profiles for rotating accretion disks if the lines originate in the very upper layers of the disk. However, Horne & Marsh (1986) noted that an anisotropy in the local emission pattern changes the global line shape: line photons are trapped in optically thick emission layers and can more easily escape in directions where there are larger Doppler gradients. The true consequence of this on the line shape will be discussed in a future paper (Hügelmeyer et al. in preparation) where we will present full 3D radiative transfer calculations in rotating accretion disks.

Figure 5: The left panel shows normalised disk ring spectra. Intensities and wavelength are offset for clarity. The right panel depicts bars corresponding in height to the contribution of each disk ring spectrum to the total spectrum. The bar for the stellar contribution has to be multiplied by a factor of . The spectra and the weights are grey-scale coded corresponding to their ring radius .
Figure 6: Temperature structures for the eight disk rings calculated for the spectral analyses. The temperature profiles marked with symbols, e. g. for disk rings with  AU and  AU, are not included in the best fit model.

3 Synthetic spectra for GQ Lup

We retrieved spectra of T Tauri stars taken with the high-resolution infrared spectrograph CRIRES at the VLT from the ESO Science Archive Facility (see Pontoppidan et al. 2008, for a description of the observations). The observations were reduced using a combination of the CRIRES pipeline and our own IDL routines. The telluric absorption lines in the spectrum were removed by using a telluric model spectrum (Seifahrt, in preparation). From the spectra taken between April 22 and April 26 2007, we selected the observation of the classical T Tauri star GQ Lup to demonstrate the applicability of our models to observations.

3.1 GQ Lup

GQ Lup is a classical T Tauri star (CTTS) which is mostly known for its recently discovered sub-stellar companion GQ Lup B (Neuhäuser et al. 2005). The activity due to ongoing accretion of the CTTS makes a well constrained determination of its physical parameters difficult. This becomes obvious regarding the differences in visual brightness of more than ( and ; Janson et al. 2006). Broeg et al. (2007) and Seperuelo Duarte et al. (2008) (BR07 and SD08 from here on) both obtained photometric and spectroscopic data to determine rotational periods and for GQ Lup A and also to derive other stellar parameters which we need as input for our models. Both authors assume an effective temperature for a K7 V star of (Kenyon & Hartmann 1995).

SD08 obtained a radius of for the K7 V star assuming a mean distance of . From evolutionary tracks they derive a mass of . Together with their rotational period of and a spectroscopically determined they find an inclination angle of .

BR07 measure a shorter photometric period of and a larger radius of from and a luminosity from evolutionary tracks. With these values and they predict an inclination angle of .

3.2 Models

We calculated small model grids for the input parameters from the publications of BR07 and SD08 (see Sec. 3.1), i. e. using , , and radii of and . The different radii lead to luminosities of and , respectively. We varied the mass accretion rate () and the Reynolds number (). For each parameter set we computed eight disk rings with . The ring radius increases exponentially to achieve a small temperature gradient from ring to ring. Figure 6 shows typical temperature structures for a set of disk rings. The decreasing temperature in the outer layers stems from the fact that our temperature is set by the gas which has a relatively low opacity at the wavelengths where stellar irradiation is the strongest. This effect becomes small at radii larger than  AU.

Since we consider a disk that has recently formed from a protostellar cloud and yet not seen significant grain growth, we assumed a standard ISM type power law grain size distribution with base size and an exponent of . The abundances follow Grevesse & Noels (1993).

In Fig. 5 we show the normalised spectra of each disk ring and how much each spectrum contributes to the total disk spectrum in the wavelength region considered. From these ring weights one can see that the outer disk radius contributes only very little to the ring integrated spectrum and that the extension to even larger disk radii would have almost no measurable influence. Furthermore, it becomes obvious that the inner rings, even though their surface area is much smaller than those of rings at larger radii, have a much higher weight. This is simply because the irradiated and viscously heated atmosphere is much warmer closer to the star and therefore has a higher flux level than that of rings further out. Furthermore, the continuum becomes optically thin for disk parts with which further decreases the flux level. The influence of the central star GQ Lup A on the total spectrum is accounted for by a blackbody spectrum of . The ratio of stellar and disk flux is depending on the input parameters for the disk. Figure 5 also shows that CO emission lines disappear at which corresponds to temperatures in the outer line emitting layers of .

3.3 Spectral fit

The comparison of our disk model spectra to the CRIRES spectrum of GQ Lup in the wavelength range indicates that the inclination of derived by SD08 is too high: the emission lines are much broader in the model than in the observation. Even if we set to large values , i. e. excluding the inner rings which have broader lines due to larger Kepler velocities, the line widths cannot be acceptably fitted. For the input parameters of BR07,

Figure 7: Spectra calculated from the same structure model but for different dust grain sizes. Neglecting dust as opacity source results in overestimation of the lines. ISM grain size distribution and a dust grain size setup with ten times smaller dust radii yield a very similar spectrum. The dust opacity is slightly reduced if we increase the grain size by a factor of ten.

we find a reasonably good fit to the CRIRES data (see Fig. 8). The line widths are best fitted with an inclination angle of , which is the lowest possible value given by the errors in BR07. The line strengths are best reproduced if we set the inner radius to and the outer radius to for a mass accretion rate of and corresponding to . The line wings of the fundamental transition lines () close to the continuum are somewhat broader in the observation than in the model. This argues for contributions of disk regions with high Kepler velocities, i. e. smaller radii . However, if we include disk rings with smaller radii, the temperature sensitive CO emission lines are too strong in the model. Furthermore, the low R- and P-branch CO lines around are overestimated by the model, or put in other words, the increment of the line strengths is too strong in the model. For our best model fit, we use a model which fits the strength of the higher order CO lines because these are more frequent in the observation accepting a less well reproduction of the low order lines.

Figure 8: CRIRES spectrum of GQ Lup (grey line) with our best fit disk model (black line). The model is calculated for a disk between and with a mass accretion rate of and a Reynolds number of . The stellar parameters are those of Broeg et al. (2007): and . For the stellar mass we adopted the value of given by Seperuelo Duarte et al. (2008). The hydrogen Pf (75) line at cannot be attributed to the disk because its Doppler width corresponds to a Kepler velocity of which lies well within our determined inner disk radius.

The part of the disk with radii has no measurable influence on the spectrum at wavelength . The strong CO fundamental transition lines () in the spectrum are well reproduced by the model. The weaker CO emission lines are also visible in the spectrum and fitted by our model. The isotope is weakly present in the model, which assumes the solar system carbon isotope ratio of 89:1, but we cannot find clear evidence for its existence in the observed spectrum. Woods & Willacy (2009) have modelled isotopic fractionation for a protoplanetary disk quite similar to our GQ Lup model, but found only mild effects on carbon monoxide, with very modest increases in the / ratio only in the outermost layers of the disk. Considering the noise in the spectrum and the blending of and lines we therefore cannot find evidence for a different isotope ratio in our observations. The maximum temperature in the surface layers of the

Figure 9: Relative abundances of the most important dust contributors to the opacity. The left plot is for the disk ring with  AU and the right one for  AU. Forsterite (MgSiO) is the most abundant dust species and strongly present in all layers. Graphite sets on around where it blocks stellar irradiation efficiently and therefore suppresses heating of deeper layers.

disk atmospheres is which is much smaller than the estimated excitation temperature of found by Najita et al. (2003) for DF Tau which has a similar spectrum but stronger emission lines.

In Fig. 7 we show the influence of the dust grain size distribution on the spectrum. All of the spectra are based on the same structure model but the emergent flux is calculated for four different dust grain size setups, i. e. no dust, ISM grain size distribution with , , and . As we can see, neglecting dust leads to a lack of opacity and too strong lines. The ISM grain size distrbution fits the data equally well as the smaller grains. The spectrum with ten times larger grains than ISM still reproduces the observation reasonably well, but lines are slightly stronger than for ISM sizes.

An interesting feature in the observed spectrum is the hydrogen Pf (75) emission line. This line has been attributed to an absorption line in the telluric standard star used for the calibration of CTTS by Najita et al. (2003). Since we used a telluric model to remove telluric absorption lines, we can relate the Pf feature unambiguously to an emission line intrinsic to GQ Lup. We measure a blue shift of for the line relative to the star and the width corresponds to a velocity of . This value is much larger than the measured by BR07. Therefore, the line cannot be attributed to the star directly. If we assume Keplerian rotation, the formal origin of the line is at and could originate from the accretion flow onto the star or a stellar or disk wind.

4 Summary & outlook

We present a new 1D structure and radiative transfer program for circumstellar disks that is an extension of the general purpose stellar atmosphere code PHOENIX. The disk is separated into concentric rings and the structure is calculated from the midplane of the disk to the top of the disk atmosphere for each ring. A combination of all rings then yields the total disk spectrum. Our program assumes a geometrically thin accretion disk geometry and a critical Reynolds number to describe the energy release in the disk due to turbulent viscosity. We have modified the hydrostatic equation, the inner and outer boundary condition of the RT equation, as well as the energy balance equation to account for the difference between star and disk atmosphere. Irradiation of the central star onto the top of the disk atmosphere is taken into account.

Our disk model spectra can reproduce high-resolution IR emission line spectra, in this case of the CTTS GQ Lup. For this particular object, we calculated a set of disk models for different stellar input parameters given by the recent publications of Broeg et al. (2007) and Seperuelo Duarte et al. (2008) and varied the mass accretion rate , the Reynolds number , and the inner and outer disk radii. We investigated the contribution of each disk ring on the total disk spectrum and showed where line and continuum radiation originates in the atmosphere. Furthermore we could show that the inclination of the warm inner disk of GQ Lup is which is much smaller than the value obtained by Seperuelo Duarte et al. (2008) and just within the uncertainties given for the inclination by Broeg et al. (2007).

Even though we could provide a reasonable model fit to the observed spectra of GQ Lup with our 1D structure and RT code, there is room for improvement. One obvious drawback of our code is the assumption of equilibrium chemistry and local thermodynamic equilibrium. This is not very likely to be valid for the cool inner and optically thin and irradiated outer disk atmospheric layers. Therefore, the consideration of departure from chemical equilibrium due to turbulent mixing or photodissociation, and radiative excitation of molecules needs to be considered in future projects. Furthermore, convective energy transport will be considered in the future in order to investigate the influence of this effect on the disk structure and spectra.

We are extending our effort to model the structure and spectra of circumstellar disks to full 3D radiative transfer calculations. This way we can relax the assumption that there is no flux exchange between matter at different radii and we will “naturally” consider the direct irradiation of the central star onto the disk and heat the very inner disk rim. The influence of the Kepler velocity field in the disk on the line profile will be investigated by means of two level atom line transfer (see Baron & Hauschildt 2007).


S.D.H. is supported by a scholarship of the DFG Graduiertenkolleg 1351 “Extrasolar Planets and their Host Stars”. A.S. acknowledges financial support from the Deutsche Forschungsgemeinschaft under DFG RE 1664/4-1.

Based on observations made with ESO Telescopes at Paranal Observatory under programme ID 179.C-0.151(A).


  • Allard et al. (2001) Allard, F., Hauschildt, P., Alexander, D., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
  • Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  • Baron & Hauschildt (2007) Baron, E. & Hauschildt, P. H. 2007, A&A, 468, 255
  • Bell & Lin (1994) Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987
  • Broeg et al. (2007) Broeg, C., Schmidt, T. O. B., Guenther, E., et al. 2007, A&A, 468, 1039
  • Cannizzo et al. (1982) Cannizzo, J. K., Ghosh, P., & Wheeler, J. C. 1982, ApJ, 260, L83
  • Carr (2007) Carr, J. S. 2007, in IAU Symposium, Vol. 243, IAU Symposium, ed. J. Bouvier & I. Appenzeller, 135–146
  • D’Alessio et al. (1998) D’Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
  • Dullemond et al. (2007) Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 555–572
  • Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  • Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  • Grevesse & Noels (1993) Grevesse, N. & Noels. 1993, in Abundances, ed. C. Jaschek & M. Jaschek (Dordrecht: Kluwer), 111
  • Hauschildt & Baron (1999) Hauschildt, P. H. & Baron, E. 1999, Journal of Computational and Applied Mathematics, 109, 41
  • Hauschildt et al. (1994) Hauschildt, P. H., Störzer, H., & Baron, E. 1994, JQSRT, 51, 875
  • Helling et al. (2008) Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854
  • Horne & Marsh (1986) Horne, K. & Marsh, T. R. 1986, MNRAS, 218, 761
  • Hubeny (1990) Hubeny, I. 1990, ApJ, 351, 632
  • Janson et al. (2006) Janson, M., Brandner, W., Henning, T., & Zinnecker, H. 2006, A&A, 453, 609
  • Kaeufl et al. (2004) Kaeufl, H.-U., Ballester, P., Biereichel, P., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. A. F. M. Moorwood & M. Iye, 1218–1227
  • Kenyon & Hartmann (1995) Kenyon, S. J. & Hartmann, L. 1995, ApJS, 101, 117
  • Klahr (2007) Klahr, H. 2007, in IAU Symposium, Vol. 239, IAU Symposium, ed. F. Kupka, I. Roxburgh, & K. Chan, 405–416
  • Kriz & Hubeny (1986) Kriz, S. & Hubeny, I. 1986, Bulletin of the Astronomical Institutes of Czechoslovakia, 37, 129
  • Lucy (1964) Lucy, L. B. 1964, SAO Special Report, 167, 93
  • Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  • Meyer & Meyer-Hofmeister (1982) Meyer, F. & Meyer-Hofmeister, E. 1982, A&A, 106, 34
  • Najita et al. (2003) Najita, J., Carr, J. S., & Mathieu, R. D. 2003, ApJ, 589, 931
  • Neuhäuser et al. (2005) Neuhäuser, R., Guenther, E. W., Wuchterl, G., et al. 2005, A&A, 435, L13
  • Pontoppidan et al. (2008) Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., et al. 2008, ApJ, 684, 1323
  • Seperuelo Duarte et al. (2008) Seperuelo Duarte, E., Alencar, S. H. P., Batalha, C., & Lopes, D. 2008, A&A, 489, 349
  • Shakura & Syunyaev (1973) Shakura, N. I. & Syunyaev, R. A. 1973, A&A, 24, 337
  • Shaviv & Wehrse (1986) Shaviv, G. & Wehrse, R. 1986, A&A, 159, L5
  • Wolf (2001) Wolf, S. 2001, PhD thesis, AA(University of Jena)
  • Woods & Willacy (2009) Woods, P. M. & Willacy, K. 2009, ApJ, accepted
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description