# Radiative transfer in circumstellar disks

###### Key Words.:

Accretion, accretion disks – Radiative transfer###### Key Words.:

Radiative transfer - Accretion, accretion disks - Methods: numerical - Techniques: spectroscopicWe 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

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

(1) |

(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

(2) |

(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

(3) |

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

(4) |

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

(5) |

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

(6) |

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

(7) |

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 .

#### 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

(8) |

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

(9) |

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

(10) |

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

(11) |

and the lower boundary condition due to symmetry conditions is

(12) |

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

(13) |

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

(14) |

and we can rewrite Eq. 13 as

(15) |

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.

(16) |

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

(17) |

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

(18) |

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

(19) |

where

(20) |

The introduction of the Eddington factors

(21) |

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

(22) | |||||

with

(23) |

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

(24) |

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

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

(26) | |||||

In Eq. 26 we assumed that the disk is axis-symmetric 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

(27) |

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.

## 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,

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.

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

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).

## 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/4-1.

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

## References

- 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