Issue 
A&A
Volume 498, Number 3, May II 2009



Page(s)  793  800  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/200811536  
Published online  19 March 2009 
Radiative transfer in circumstellar disks
I. 1D models for GQ Lupi
S. D. Hügelmeyer^{1}  S. Dreizler^{1}  P. H. Hauschildt^{2}  A. Seifahrt^{1}  D. Homeier^{1}  T. Barman^{3}
1  Institut für Astrophysik, GeorgAugustUniversität
Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
2 
Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany
3 
Lowell Observatory, 1400 W Mars Hill Rd, Flagstaff, AZ 86001, USA
Received 17 December 2008 / Accepted 10 March 2009
Abstract
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 highresolution VLT CRIRES
spectra of the T Tauri star GQ Lup.
Key words: radiative transfer  accretion, accretion disks  methods: numerical  techniques: spectroscopic
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 10years 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 LyndenBell & Pringle (1974) or models using the diffusion approximation (e.g. Cannizzo et al. 1982; Meyer & MeyerHofmeister 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 midinfrared. Highresolution 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 shortperiod 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 Sect. 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 Sect. 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 welltested model atmosphere packagePHOENIX
(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
phaseequilibrium 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 S=1; 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
wellmixed disk atmosphere. In the lowdensity outer layers,
nonthermal 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.
Figure 1: Disk ring structure as adopted for our calculations. The radius of the rings increases exponentially. The left panel shows a faceon view of a disk, the right one a vertical cut viewed edgeon (height is not to scale). The dotted lines are the radii R 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. 

Open with DEXTER 
In our models, the disk region considered is divided into rings (see Fig. 1) and a planeparallel 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 H is much smaller than the disk ring radius
R. This assumption decouples the treatment of vertical and radial
disk structure because the vertical structure is in quasistatic
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 T<200 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 dissipation and angular momentum transport but its effect on the thermal structure
of the disk cannot be easily described or parametrised. Even though
temperatures T>1000 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
often modelled as an ``alphaviscosity'' resulting in a
verticallyaveraged viscosity
(Shakura & Syunyaev 1973) where is the angular momentum transfer efficiency, the sound speed, and H the pressure scale height. Alternatively, a turbulent viscosity can be modelled assuming a critical Reynolds number Re
(LyndenBell & 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 newDISK
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 NL is the number of layers, which
will be kept fixed also for later calculations with the same input
parameters. The m scale is equally spaced on a logarithmic grid
between the column mass at the midplane of the disk
and the outer value m_{1}, which is an input parameter. This spacing is done for points, while for the remaining layers each point j, , is positioned halfway between M_{0} and m_{j1}, with . This is done to provide numerical stability in the iterative process because of the steep mgradient near the midplane. A typical value for is 6 or 12 for NL=64.
In a next step, we calculate the depthdependent 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. LyndenBell & 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 m into a height above the midplane of the disk z.
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 Tstructure is consistent for each layer. In Eq. (6) we set equal to 1 and is the optical thickness at the midplane of the disk (see Sect. IIIb 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 iswhere 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 F are of the order of 10^{2}.
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öldLucy temperature correction is reached after iterations. 

Open with DEXTER 
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 selfgravitation 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 diskwhere P 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 z (which is not known apriori) and introduce the sound speed , which is taken from the previous iteration and only depends on the temperature T. Therefore, the error in P and cancel out and we obtain a more stable version of Eq. (8). This secondorder differential equation is then decomposed into two coupled firstorder differential equations. With the inner boundary condition (symmetry at the midplane) and the outer boundary condition P(m_{1}) 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 planeparallel radiative transfer equationfor a given number of Gaussian angles , , with corresponding quadrature weights a_{i} from the midplane (z=0, m=M_{0}) to the top of the disk ( , m=m_{1}). 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 S is the source function and is the thermal coupling parameter. The operator is split according to
and we can rewrite Eq. (13) as
The nonlocal 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öldLucy class temperature correction scheme (Lucy 1964) by integrating the second frequencyintegrated moment of the radiative transfer equation
where
The introduction of the Eddington factors
and insertion of the first frequencyintegrated moment of the radiative transfer equation gives the temperature correction law
with
the absorption mean and Planck mean opacity. An example iteration history is shown in Fig. 2.
Figure 3: Irradiation geometry as adopted for our calculations. We consider a star with radius at a distance R from a disk ring with a height and a slope of the disk surface . For each Gaussian angle and its corresponding integration weight a_{i} 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 . 

Open with DEXTER 
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 R by computing the height of the disk according to our start model (see Sect. 2.1) at with . For each Gaussian angle and corresponding integration weight a_{i} the projected surface fraction on the star A_{i} 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 towhere NS is the number of surface segments for Gauss angle , r_{j} is the distance from the area element A_{j} 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.
Figure 4: Optical depth structure for a disk ring model with , , , , and . The left plot shows the run of with height (z/R=0 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 temperaturediagram showing at which temperatures and column masses the optical depth for the line (dashed line) and the continuum (dotted line) become unity. 

Open with DEXTER 
2.2.5 Spectrum
Having determined a set of selfconsistent disk ring structures, we calculate a last iteration with a larger number of Gaussian angles (usually NG=16) and a fine wavelength spacing (up to in the wavelength range of interest) to obtain a highresolution model spectrum. Since the Keplerian rotation velocity and the inclination angle i of the disk are not considered in the single ring spectra yet, we compute a combined spectrum according toIn Eq. (26) we assumed that the disk is axissymmetric and that the intensity is constant for all radii between inner and outer 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 iand for a set of disk ring segments with azimuthal angle . We use NA=100 disk ring segments, i.e. 100 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.
3 Synthetic spectra for GQ Lup
We retrieved spectra of T Tauri stars taken with the highresolution 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 substellar 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 L from evolutionary tracks. With these values and they predict an inclination angle of .
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 25. The spectra and the weights are greyscale coded corresponding to their ring radius R. 

Open with DEXTER 
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 R=0.031 AU and R=0.045 AU, are not included in the best fit model. 

Open with DEXTER 
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. 

Open with DEXTER 
3.2 Models
We calculated small model grids for the input parameters from the publications of BR07 and SD08 (see Sect. 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 R 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 R=0.065 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 3.5. 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 .
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. 

Open with DEXTER 
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, 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 (v=10) 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 R. However, if we include disk rings with smaller radii, the temperature sensitive v=21 CO emission lines are too strong in the model. Furthermore, the low J R and Pbranch CO v=10 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 v=10 lines because these are more frequent in the observation accepting a less well reproduction of the low order lines.The part of the disk with radii has no measurable influence on the spectrum at wavelength . The strong CO fundamental transition lines (v=10) in the spectrum are well reproduced by the model. The weaker v=21 CO emission lines are also visible in the spectrum and fitted by our model. The ^{13}CO 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 ^{12}CO/ ^{13}CO ratio only in the outermost layers of the disk. Considering the noise in the spectrum and the blending of ^{12}CO and ^{13}CO lines we therefore cannot find evidence for a different isotope ratio in our observations. The maximum temperature in the surface layers of the 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 v=21 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 Pffeature 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.
Figure 9: Relative abundances of the most important dust contributors to the opacity. The left plot is for the disk ring with R=0.094 AU and the right one for R=0.290 AU. Forsterite (Mg_{2}SiO_{4}) 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. 

Open with DEXTER 
4 Summary and 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 highresolution 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 Re, 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 R 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).
Acknowledgements
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/41. Based on observations made with ESO Telescopes at Paranal Observatory under programme ID 179.C0.151(A).
References
 Allard, F., Hauschildt, P., Alexander, D., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357 [NASA ADS] [CrossRef] (In the text)
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] (In the text)
 Baron, E., & Hauschildt, P. H. 2007, A&A, 468, 255 [NASA ADS] [CrossRef] [EDP Sciences]
 Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] (In the text)
 Broeg, C., Schmidt, T. O. B., Guenther, E., et al. 2007, A&A, 468, 1039 [NASA ADS] [CrossRef] [EDP Sciences]
 Cannizzo, J. K., Ghosh, P., & Wheeler, J. C. 1982, ApJ, 260, L83 [NASA ADS] [CrossRef]
 Carr, J. S. 2007, in IAU Symp. 243, ed. J. Bouvier, & I. Appenzeller, 135 (In the text)
 D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [NASA ADS] [CrossRef] (In the text)
 Dullemond, C. P., Hollenbach, D., Kamp, I., & D'Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 555 (In the text)
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] (In the text)
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
 Grevesse, N., & Noels. 1993, in Abundances, ed. C. Jaschek, & M. Jaschek (Dordrecht: Kluwer), 111 (In the text)
 Hauschildt, P. H., & Baron, E. 1999, J. Comput. Appl. Math., 109, 41 [NASA ADS] [CrossRef] (In the text)
 Hauschildt, P. H., Störzer, H., & Baron, E. 1994, JQSRT, 51, 875 [NASA ADS] (In the text)
 Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854 [NASA ADS] [CrossRef] (In the text)
 Horne, K., & Marsh, T. R. 1986, MNRAS, 218, 761 [NASA ADS] (In the text)
 Hubeny, I. 1990, ApJ, 351, 632 [NASA ADS] [CrossRef] (In the text)
 Janson, M., Brandner, W., Henning, T., & Zinnecker, H. 2006, A&A, 453, 609 [NASA ADS] [CrossRef] [EDP Sciences]
 Kaeufl, H.U., Ballester, P., Biereichel, P., et al. 2004, in SPIE Conf. Ser. 5492, ed. A. F. M. Moorwood, & M. Iye, 1218 (In the text)
 Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117 [NASA ADS] [CrossRef] (In the text)
 Klahr, H. 2007, in IAU Symp., 239, ed. F. Kupka, I. Roxburgh, & K. Chan, 405 (In the text)
 Kriz, S., & Hubeny, I. 1986, Bull. Astron. Inst. Czechoslovakia, 37, 129 [NASA ADS] (In the text)
 Lucy, L. B. 1964, SAO Special Report, 167, 93 [NASA ADS] (In the text)
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] (In the text)
 Meyer, F., & MeyerHofmeister, E. 1982, A&A, 106, 34 [NASA ADS]
 Najita, J., Carr, J. S., & Mathieu, R. D. 2003, ApJ, 589, 931 [NASA ADS] [CrossRef] (In the text)
 Neuhäuser, R., Guenther, E. W., Wuchterl, G., et al. 2005, A&A, 435, L13 [NASA ADS] [CrossRef] [EDP Sciences]
 Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., et al. 2008, ApJ, 684, 1323 [NASA ADS] [CrossRef] (In the text)
 Seperuelo Duarte, E., Alencar, S. H. P., Batalha, C., & Lopes, D. 2008, A&A, 489, 349 [NASA ADS] [CrossRef] [EDP Sciences]
 Shakura, N. I., & Syunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS]
 Shaviv, G., & Wehrse, R. 1986, A&A, 159, L5 [NASA ADS]
 Wolf, S. 2001, Ph.D. Thesis, AA, University of Jena (In the text)
 Woods, P. M., & Willacy, K. 2009, ApJ, accepted (In the text)
All Figures
Figure 1: Disk ring structure as adopted for our calculations. The radius of the rings increases exponentially. The left panel shows a faceon view of a disk, the right one a vertical cut viewed edgeon (height is not to scale). The dotted lines are the radii R 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. 

Open with DEXTER  
In the text 
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öldLucy temperature correction is reached after iterations. 

Open with DEXTER  
In the text 
Figure 3: Irradiation geometry as adopted for our calculations. We consider a star with radius at a distance R from a disk ring with a height and a slope of the disk surface . For each Gaussian angle and its corresponding integration weight a_{i} 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 . 

Open with DEXTER  
In the text 
Figure 4: Optical depth structure for a disk ring model with , , , , and . The left plot shows the run of with height (z/R=0 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 temperaturediagram showing at which temperatures and column masses the optical depth for the line (dashed line) and the continuum (dotted line) become unity. 

Open with DEXTER  
In the text 
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 25. The spectra and the weights are greyscale coded corresponding to their ring radius R. 

Open with DEXTER  
In the text 
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 R=0.031 AU and R=0.045 AU, are not included in the best fit model. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
Figure 9: Relative abundances of the most important dust contributors to the opacity. The left plot is for the disk ring with R=0.094 AU and the right one for R=0.290 AU. Forsterite (Mg_{2}SiO_{4}) 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. 

Open with DEXTER  
In the text 
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.