© ESO, 2010

1. Introduction

According to the general consensus, planets and planetary systems form in circumstellar disks. These disks consist of gaseous and solid dust material. There are two major mechanisms proposed for the formation of giant planets. The first is the disk instability, which requires a gravitationally unstable disk in which giant planets form by direct collapse of the gas as a consequence of its self-gravity, originally suggested by Kuiper (1951); Cameron (1978) and more recently by Boss (2001). Later on, the gas giant may collect dust, which after settling to its center, forms a solid core. The second mechanism is the core-accretion process (Bodenheimer & Pollack 1986; Pollack et al. 1996), which is the final stage of planet formation in the planetesimal hypothesis (Safronov 1972). In this hypothesis dust coagulates first and forms planetesimals (meter to kilometer-sized objects). The subsequent growth of planetesimals consists of two stages, the runaway (Wetherill & Stewart 1989) and the oligarchic growth (Kokubo & Ida 1998), in which the formed bodies are massive enough to increase their masses by gravity-assisted collisional accretion processes. In this way planetary embryos and, through their consecutive collisions, protoplanets, terrestrial planets, and planetary cores of giant planets can be formed. During the gas accretion phase the precursor of the proto giant planet reaches a critical mass (where the planetary envelope contains more material than its core), and at about the order of 10   M_⊕ starts to contract, causing an increased accretion rate, which in turn raises radiative energy losses resulting in a runaway gas accretion. In this model the planetary core rapidly builds up a massive envelope of gas from its surrounding disk in less than 10   Myr, and eventually a gas giant forms. Because the disk dispersal time is presumably about 5   Myr (Haisch et al. 2001; Hillenbrand 2005), this model should also incorporate planetary migration to explain giant planets observed well inside the snow line, see Alibert et al. (2004).

Both planetary formation processes have their weak points. Disk instability assumes a very massive disk that is gravitationally unstable or only marginally stable. Another as yet unsolved problem is that during the gravitational collapse an effective cooling mechanism should work to ensure the formation of gravitationally bound clumps. While the radiative cooling is not effective enough and the recently invoked thermal convection is also ineffective, the disk instability as a common planet-formation process is questionable (Klahr 2008).

Regarding the giant planet formation in the planetesimal hypothesis, Mordasini et al. (2009) convincingly presented the success of the core accretion model to produce giant planets within the disk lifetime in a planet population synthesis model, but the artificial slow-down of the type I migration was necessary. Type I migration has such a short time scale (Ward 1997) that the planet inevitably will be engulfed by the host star within the disk lifetime. Another unsolved issue of the planetesimal hypothesis is the so-called “meter-sized barrier problem”: i.e. the formation of meter-sized bodies is impeded by their quick inward drift to the star (Weidenschilling 1977), and by mutual disruptive collisions due to their high relative velocities above 1m/s (Blum & Wurm 2008). If no other physical processes are taking place, these two mechanisms would impede the formation of the meter-sized and consequently the larger planetesimals.

A wealth of information about the planet formation process would clearly be provided if newly formed planets, or at least their signatures, could be observed in still gas-rich protoplanetary disks. There are already attempts to discover signatures of giant planet-formation assuming the disk instability mechanism in the works of Narayanan et al. (2006) and Jang-Condell & Boss (2007). In the first study, for instance, the authors suggest observations of HCO⁺ emission lines as a tracer of the accumulation of large clumps of cold material, arguing that with high-resolution millimeter and sub-millimeter interferometers the cold clumps can be directly imaged. The observability of an already formed planet embedded in a circumstellar disk has also been studied recently by Wolf & D’Angelo (2005) and Wolf et al. (2007). In these works the authors investigated whether the influence of a forming planet on the protoplanetary disk can be detected by analyzing the spectral energy distribution (SED) coming from the star-disk system. It was found that by studying only the SED it is impossible to infer the presence of an embedded planet. On the other hand, it was also concluded that the dust re-emission in the hot regions near the gap could certainly be detected and mapped by ALMA for the nearby (less than 100   pc in distance) and approximately face-on protoplanetary disks.

Another interesting attempt was presented by Clarke & Armitage (2003), where the authors have investigated the possibility of detection of giant planets by excess CO overtone emission originating from the planetary accretion inflow in FU Ori type objects. Nevertheless, the periodic line profile distortions are observable only in tight systems in which the giant planet is orbiting within 0.25   AU. Another possibility to detect embedded giant planets is based on the radial velocity measurements, which favor high-mass, short-period companions. The detectability limit of an embedded giant planet by radial velocity measurements is  ~10   M_J or  ~2   M_J (M_J is the mass of Jupiter) orbiting at 1 AU or 0.5 AU, respectively, due to heavy variability in the optical spectra of young (<10   Myr) stars (e.g., Paulson & Yelda 2006; Prato et al. 2008). Note that although there were attempts to detect planets by radial velocity measurements in young systems, no firm detection has yet been repeated see, e.g., Setiawan et al. (2008) for TW Hya b, debated later by Huélamo et al. (2008). The discovery of planetary systems by direct imaging nowadays is restricted to large separations  >10−100   AU (Veras et al. 2009) and favors higher mass planets. As of today nine exoplanetary systems have been directly imaged in the optical or near-infrared band, but only five of them are still embedded: GQ Lup (Neuhäuser et al. 2005) and CT Cha (Schmidt et al. 2008), but they are rather hosting a stellar companion as the companion masses are, because high as 21.5   M_J and 17   M_J, respectively; 2M1207 (Chauvin et al. 2004, 2005a) hosting a 4   M_J mass planet; UScoCTIO 108 (Kashyap et al. 2008) hosting a 14   M_J mass planet; βPic (Lagrange et al. 2009a,b) hosting an 8   M_J mass planet. AB Pic (Chauvin et al. 2005b), HR 8799 (Marois et al. 2008) (triple system), SR 1845 (Biller et al. 2006) and Fomalhaut (Kalas et al. 2008) are mature planetary systems with ages of about 30 Myr, 60 Myr, 100 Myr, and 200 Myr. Note that centimeter wavelength radio observations of a very young (<0.1   Myr) disk around HL Tau, presented by Greaves et al. (2008), revealed the possibility of a 12 Jupiter mass giant planet being in formation at  ~75   AU.

We investigate a new possibility to detect young giant planets embedded in protoplanetary disks of T Tauri stars via line profile distortions in its near-IR spectra. Contrary to the studies mentioned above, instead of dealing with the disk’s density perturbations, we investigate the gas dynamics in the disk, particularly near the gap opened by the planet. We have found that there is a reasonable possibility to discover giant planets embedded in the protoplanetary disk of young stars with our approach.

Our idea has been inspired by recent findings of Kley & Dirksen (2006), who showed that a massive giant planet embedded in a protoplanetary disk produces a transition of the disk from the nearly circular state into an eccentric state. The mass limit for this transition is 3   M_J for a disk with ν = 10^-5 uniform kinematic viscosity (measured in dimensionless units, see details in Sect. 3.1). The most visible manifestation of the disk’s eccentric state is that the outer rim of the gap opened by the giant planet has an elliptic shape. After performing a series of hydrodynamical simulations with the code FARGO (Masset 2000), we can also confirm the results of Kley & Dirksen (2006) even using α-type viscosity. Note that beside the elliptic geometry that is clearly visible in the disk surface-density distribution, the giant planet also disturbs the motion of the gas parcels near the gap. Using CO as tracer molecule of gas, which is excited to  ~3−5   AU (depending on the stellar mass) in the disk atmosphere (Najita et al. 2007), we investigate how the CO emission line profiles at 4.7   μm are distorted by the gas dynamics. With our semi-analytical synthetic spectral model we explore the influence of the mass of the giant planet, the disk inclination angle, the mass of the hosting star, the orbital distance of planet, and finally, the size of inner cavity on the line profile distortions.

The paper is structured as follows: in Sect. 2 we review the basic physics of our synthetic spectral line models. In Sect. 3 we present disk models and details of numerical simulations. Our results of the CO ro-vibrational line profile distortions, the detectability of an embedded giant planet, and its observability constraints are presented in Sect. 4. The paper closes with the discussion of the results and concluding remarks.

2. Synthetic spectral line model

To calculate the CO ro-vibrational spectra emerging from protoplanetary disks we developed a semi-analytical line spectral model. The thermodynamical model of the disk is based on the two-layer flared disk model, originally proposed by Chiang & Goldreich (1997), which describes the heating by stellar irradiation in a disk with a flared geometry. In our model we use the flared disk approximation, and assume the disk to consist of two layers: an optically thick interior, producing continuum radiation and an optically thin atmosphere, producing line emission or absorption in the spectra¹. The boundary between the two layers lies where the dust becomes optically thick in the visual wavelengths, i.e. at optical depth τ_V = 1 along the line of the incident stellar irradiation. The irradiation from the boundary layer formed at the inner edge of the disk (Popham et al. 1993) is neglected though because its existence had not yet been confirmed conclusively by observations yet. Nevertheless, the emission from the inner disk rim is substantial, resulting in strong near-IR bump in Herbig Ae/Be (Natta et al. 2001) and weaker one in T Tauri SEDs (Muzerolle et al. 2003). In addition it was found by Monnier et al. (2006) that the simple vertical wall assumption for the inner rim is inconsistent with the interferometric measurements. Indeed, the shape of the disk inner rim is presumably rounded-off, resulting in less inclination-angle-dependent near-IR excess emission (Isella & Natta 2005). A schematic representation of the disk model is shown in Fig. 1. If the disk atmosphere is superheated by the stellar irradiation with respect to its interior, we may expect emission lines. Conversely, if the disk interior is heated above the temperature of the disk atmosphere by viscous dissipation, for example due to an abrupt increase in accretion rate, spectral lines are expected in absorption.

The vertical structure of a geometrically thin disk can be derived by considering vertical hydrostatic equilibrium (Shakura & Sunyaev 1973). Ignoring any contribution from the gravitational force of the disk itself, the vertical density profile would be set by the equilibrium of gas pressure, the stellar gravitation and centrifugal force owing to orbital motion, resulting in a Gaussian vertical density profile. Because the thickness of the disk atmosphere in the double layer model is principally determined by the grazing angle of the stellar irradiation, which is narrow (δ(R) ≪ 1) in our computational domain (R ≤ 5   AU), we assume that the disk atmosphere has uniform vertical density distribution. Note that the surface density in the disk atmosphere according to Eq. (D.2) does not depend on the surface density (Σ(R)) distribution by definition.

The expected flux emerging from the double-layer disk is the result of the radiation of gas emission from the optically thin disk atmosphere above the continuum emission of dust in the optically thick disk interior and inner rim. Because the stellar continuum contributes considerably to the line flux in the near-IR band compared to the disk continuum, especially in disks with an inner cavity, we have to also take into account the stellar continuum in the total flux. In our model the stellar radiation is taken to be blackbody radiation of the stellar effective surface temperature. Because the gas temperature is regulated by collisions with dust grains and stellar X-ray heating (which is significant for a young star), the gas and dust components may be thermally uncoupled in the tenuous disk atmosphere below a critical density of n_cr ~ 10¹³−10¹⁴   cm^-3 (Chiang & Goldreich 1997). As a result the gas temperature can be as high as  ~4000−5000   K in a region where the column density is below  ~10²¹   cm^-2 (Glassgold et al. 2004). Because the gas column density is  ~10²²   cm^-2 in our model disk atmosphere, we adopt thermal coupling of dust and gas, i.e. T_g(R) = T_d(R) = T_atm,irr(R), where T_atm,irr(R) is the atmospheric temperature defined by Eq. (A.5). Although Glassgold et al. (2004) and Kamp & Dullemond (2004) suggest the existence of an overheated layer above the superheated disk atmosphere, we did not consider this for simplicity’s sake. In this model, assuming local thermodynamic equilibrium, the monochromatic intensity at frequency ν emitted by gas parcels at a distance R to the star and an azimuthal angle φ observed at inclination angle i (i = 0° meaning face on disk) can be given by $\begin{array}{ccc}\mathit{I}\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}& \mathrm{=}& \mathit{B}\mathrm{\left(}\mathit{\nu ,}{\mathit{T}}_{\mathrm{int}}\mathrm{\right(}\mathit{R}\mathrm{\left)}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}}\\ & \mathrm{+}& \mathit{B}\mathrm{\left(}\mathit{\nu ,}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right(}\mathit{R}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}}\mathrm{)}\mathit{,}\end{array}$(1)where τ(ν,R,φ,i) is the monochromatic optical depth of the disk atmosphere at a frequency ν along the line of sight, B(ν,T_int(R)) and B(ν,T_atm(R)) are the Planck functions of the dust temperature T_int(R) in the disk interior, and the gas temperature T_atm(R) in disk atmosphere, respectively. Regarding the disk inner edge, which has a temperature T_rim(R₀), it is assumed that it radiates as a blackbody (I_rim(ν) = B(ν,T_rim(R₀))), see Appendix B for details. The total flux emerging from the protoplanetary disk seen by inclination angle i at frequency ν can thus be given as $\begin{array}{ccc}{\mathit{F}}_{\mathrm{disk}}\mathrm{\left(}\mathit{\nu ,i}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}{\mathit{R}}_{\mathrm{1}}\\ {\mathit{R}}_{\mathrm{0}}\end{array}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{I}\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}\frac{\mathit{RdRd\phi }}{{\mathit{D}}^{\mathrm{2}}}\cos \mathrm{\left(}\mathit{i}\mathrm{\right)}\\ & \mathrm{+}& \frac{{\mathit{I}}_{\mathrm{rim}}\mathrm{\left(}\mathit{\nu ,R}\mathrm{\right)}}{{\mathit{D}}^{\mathrm{2}}}{\mathit{A}}_{\mathrm{rim}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\cos \mathrm{(}\mathrm{90}\mathrm{-}\mathit{i}\mathrm{)}\mathit{,}\end{array}$(2)where R₀ and R₁ are the disk inner and outer radii and D is the distance of the source to the observer. In Eq. (2) A_rimcos(90 − i) is the visible area of the disk inner rim. If the disk inner rim were a perfect vertical wall, its surface area would be , where h(R₀) is the disk aspect ratio at the inner edge, which is taken to be 0.05. In order to be consistent with the observations regarding the inclination-angle dependence of the rim flux (Isella & Natta 2005), the rim is taken to be a vertical wall seen under a constant 60 degree inclination angle.

To calculate the synthetic spectra numerically, we first set up the two-dimensional computational domain with N_R logarithmically distributed radial and N_φ equidistantly distributed azimuthal grid cells. For a given set of stellar parameters, shown in Table 1, first the grazing angle of the incident stellar irradiation is determined according to Eq. (A.2). The temperature distributions in the disk atmosphere and interior are calculated according to Eq. (A.5) and Eqs. (A.7, A.9, B.2, B.3), respectively. After determining the dust surface-density of the disk atmosphere with Eq. (D.2), the optical depth along the line of sight is calculated by Eq. (D.4) using the gas opacity given by Eq. (E.1) and considering Eqs. (E.4)–(E.8) to calculate appropriate intrinsic line profiles and Doppler shifts. In order to determine the mass absorption coefficient of the molecule we used the transition data provided by Goorvitch (1994), such as transition probability A, lower and upper state energy E_l,E_u and statistical weight g_u. Knowing the optical depth (τ(ν,R,φ,i)), the temperature in disk atmosphere and interior in the 2D computational domain the expected flux along the line of sight can be calculated by applying Eqs. (1) and (2) at each frequency in the vicinity of the investigated transition.

3. Disk models

In this section we describe in detail our disk models that fed into the synthetic spectral line model. As a first step, we calculated the disk surface-density and velocity distribution perturbed by a massive embedded planet with the publicly available hydrodynamical code FARGO of Masset (2000). A very important difference to the planet-free case is that the disk surface-density distribution is very heavily modified. It shows a clear elliptic character of the gap opened by the giant planet. According to Kley & Dirksen (2006), D’Angelo et al. (2006), and also to our calculations (see below), the disk becomes eccentric after several hundred orbits of the embedded planet. With regard to our synthetic spectral line model, the most relevant are those effects that are originating from the gas dynamics caused by the giant planet. It is well known that a pure Keplerian motion of the gas parcels on circular orbits results in symmetric double-peaked emission line profiles (Horne & Marsh 1986). As the CO emission is strongly depressed 2−3   AU in our models, the origin of the double-peaked profiles is clear. Any deviation of the gas dynamics from the pure rotation will break this symmetry, resulting in the distortion of line profiles. If the distorting effect is strong enough, the line profile distortion could be significant to be detected by high-resolution spectroscopic observations.

Indeed several T Tauri stars show symmetric but centrally peaked CO line profiles that can be explained by radially extent CO emission overriding the double peaks (Najita et al. 2003; Brittain et al. 2009). Stronger CO emission can be expected if UV fluorescence plays a role (Krotkov et al. 1980), or the gas temperature is not well coupled to the dust (Kamp & Dullemond 2004; Glassgold et al. 2004). In the tenuous region above the disk atmosphere the gas can be significantly hotter than is predicted by the double-layer model. As a consequence the slowly rotating distant disk parcels could produce a substantial contribution to the low-velocity part of the line profile, resulting in a centrally peaked profile. Nevertheless, for simplicity’s sake here we do not consider any of the above mentioned effects.

3.1. Hydrodynamical setup

We adopt dimensionless units in hydrodynamical simulations, for which the unit of length and mass is taken to be the orbital distance of the planet, and the mass of the central star, respectively. The unit of time, t₀, is taken to be the reciprocal of the orbital frequency of the planet, resulting in t₀ = 1/2π, setting the gravitational constant to unity. In each disk model we assigned four planetary masses to the embedded planet, which are in increasing order q = m_pl/m_∗ = 0.001, 0.003, 0.005, and 0.008, where m_pl and m_∗ are the planetary and stellar masses. Regarding the disks geometry, the flat thin disk approximation was assumed with an aspect ratio of H(R)/R = 0.05. The disk extends between 0.2−5 dimensionless units. This computational domain was covered by 256 radial and 500 azimuthal grid cells. The radial spacing was logarithmic, while the azimuthal spacing was equidistant. This results in practically quadratic grid cells, because the approximation ΔR ~ RΔφ is valid at each radius. The disks are driven by α-type viscosity (Shakura & Sunyaev 1973), which is consistent with our thermal disk model with an intermediate value of α, which was taken to be 1 × 10^-3. Note that assuming a cold disk (therefore virtually non-ionized) the viscosity generated by magnetohydrodynamical turbulence is hard to explain, although, according to Mukhopadhyay (2008), the transient growth of two or three-dimensional pure hydrodynamic elliptic-type perturbations could result in α ≃ 10^-1−10^-5, depending on the disk scale height (decreasing α with increasing scale height). The disk’s initial surface density profile is given by a power law Σ(R) = Σ₀R^−1/2, where the surface density at 1 distance unit Σ₀ is taken to be 2.15 × 10^-5 in dimensionless units. For simplicity, a locally isothermal equation of state is applied for the gas, and the disk self-gravity is neglected. While the inner boundary of the disk is taken to be open, allowing the disk material to leave the disk on accretion timescale, the outer boundary is closed, i.e. no mass supply is allowed.

Below we assign physical units to the dimensionless quantities. If the stellar mass is assumed to be m_∗ = 1   M_⊙, the planetary masses in our models correspond to 1   M_J, 3   M_J, 5   M_J, and 8   M_J. Thanks to the dimensionless calculations it is possible to scale the results to different stellar masses. For different stellar masses, the mass of the giant planet scales with the stellar mass, i.e. m_pl = 1 × 10^-3   m_∗, 3 × 10^-3   m_∗, 5 × 10^-3   m_∗, 8 × 10^-3   m_∗, and the stellar masses are taken to be m_∗ = 0.5   M_⊙, 1   M_⊙, 1.5   M_⊙. Furthermore, the distance unit is set to the distance of the planet to the star, i.e., the planet orbits at 1 AU. Below we extended our models to tight and wide systems, where the giant planets are orbiting at 0.5 and 2 AU, respectively. In 1   M_⊙ stellar mass disk models the surface density at 1 AU is Σ₀ = 2.15 × 10^-5   M_⊙AU^-2 ≃ 191   g/cm² at the beginning of simulation, which is also scaled with the mass of the central star. But note that the solution to the hydrodynamical system of equations is independent of Σ₀, if there is no back-reaction to the planet, see the vertically integrated Navier-Stokes equations in Kley (1999). For 0.5   M_⊙, 1   M_⊙ and 1.5   M_⊙ stellar mass model the disk mass in the computational domain (0.2   AU ≤ R ≤ 5   AU) is 5 × 10^-4   M_⊙, 1 × 10^-3   M_⊙ and 1.5 × 10^-3   M_⊙. These values refer to the mass of the inner disk only, as the whole disk may extend to several tens or hundreds of AU, and the total disk mass is in the conventional range of 0.01−0.1   M_⊙. More details about the models can be found in Table 1.

3.2. Synthetic spectral calculation setup

We calculate the spectral line of a fundamental band CO transition (V = 1 → 0P(10), λ₀ = 4.7545   μm), which is not blended by the higher excitation V = 2 → 1 and V = 3 → 2 lines, although in our thermal model the temperature is not high enough to excite the higher vibrational levels at all.² The CO rotational and vibrational level populations were calculated in local thermodynamical equilibrium. The line profiles were calculated in 200 wavelengths and in the same numerical domain that was used in hydrodynamical simulations. According to our preliminary tests on “circularly Keplerian”³ disks, the CO fundamental band spectra did not show significant changes with increasing outer disk boundary. This can be explained by the fact that a substantial part of the CO fundamental band emission is arising in the inner disk, consequently the outer part of the disk (R > 3−5   AU, depending on the stellar mass) does not contribute to the disk emission at 4.7 micron. Contrary to this, the distance of the inner boundary of the disk to the central star heavily influences the CO fundamental band spectra, i.e. the closer the disk inner boundary to the star the line-over-continuum is weaker and broader. This is because the disk innermost regions give stronger contribution to the continuum, which in turn weakens the CO emission normalized to continuum, although the CO emission itself is also getting stronger due to an increased amount of hot CO. The line broadening is a natural consequence of the larger orbital velocity of gas orbiting closer to the star.

Because we investigate disks with already formed planets, it is reasonable to assume that a cavity is formed at the very inner disk like the ones found in disks of several T Tauri stars (Akeson et al. 2005; Eisner et al. 2005, 2007; Salyk et al. 2009). If dust evaporation is takes place close to the star, the CO will be depleted too, due to the disappearance of the dust which protects the CO molecules against UV photo dissociation. In this way, the observed inner radius of CO (R_CO) should intuitively be similar or slightly larger than the dust sublimation radius R_sub. Indeed, R_CO was found to be larger than R_sub by a factor of a few in several disks, e.g., see Fig. 5. of Eisner et al. (2007) or Fig. 10 in Salyk et al. (2009). Contrary to this, Carr (2007) found that the CO inner radius is  ~0.7 that of the dust. If the gas disk extends as far as 0.05 AU, there would be an additional broad component to the emission that is not modeled here. Taking these argumentations into account it is reasonable to assume that the CO inner radius is set by the dust sublimation. The dust sublimation radius R_sub is given by ${\mathit{R}}_{\mathrm{sub}}\mathrm{\simeq }\mathrm{0.03}{\left(\frac{{\mathit{T}}_{\mathrm{\ast }}}{{\mathit{T}}_{\mathrm{sub}}}\right)}^{\mathrm{5}\mathit{/}\mathrm{2}}\frac{{\mathit{R}}_{\mathrm{\ast }}}{{\mathit{R}}_{\mathrm{\odot }}}\mathrm{AU}\mathit{,}$(3)where Eq. (A.5) with standard dust opacity law, β = 1 (Rodmann et al. 2006) was used. Assuming T_sub = 1500   K for dust sublimation temperature, T_∗ ~ 4000   K and R_∗ = 1.8   R_⊙ for stellar surface effective temperature and radius, appropriate for a 2.5 Myr solar mass star, we obtain R_sub = 0.05   AU. To keep this simple, we set the disk inner boundary (namely the CO inner radius) fixed to R_CO ≡ 4   R_sub = 0.2   AU, which is an acceptable value in disks hosted by T Tauri stars with 3500–5000 K surface temperature and 1.4−2.2   R_⊙ radius. Note that R_sub was fixed throughout our models (except the ones where we investigated the effect of the size of the inner cavity) to make them comparable, neglecting the effect of the change in stellar luminosity on the R_sub.

It is reasonable to assume that a disk that is 2.5   Myr old contains a considerable amount of gas and tracer CO also in the inner disk, hence the stellar age is taken to be 2.5   Myr in all models. Note that according to Haisch et al. (2001), half of the observed young stars in nearby embedded clusters (NGC 2024, Trapezium, IC 348, NGC 2264, NGC 2362, NGC 1960), with ages of about 3   Myr, show near-IR excess, which is related to the excess emission of dust in disks. Furthermore, Hillenbrand (2005) also concluded that the median lifetime of the optically thick inner disk is between 2−3   Myr.

In order to simplify the models it is assumed that the dust consists of pure silicates with 0.1   μm grain size (Draine & Lee 1984). According to this the mass absorption coefficient of the dust is taken to be 2320   cm²/g at visual wavelength. Note that neither the effect of the coagulation nor the settling out of tenuous atmosphere of grains with a large size are taken into account. The size distribution and therefore the overall mass absorption coefficient of the dust is taken to be the same in the interior and the atmosphere. In this way the dust opacity is slightly overestimated. According to Eq. (D.4) the atmospheric gas density and thus gas line emission strength is slightly underestimated. The dust-to-gas and CO-to-gas mass ratios are assumed to be constant throughout the disk and are taken to be 10^-2 and 4 × 10^-4 (measured per gram of gas+dust mass) in the disk atmosphere and interior, respectively.

For the planet-free disks circular Keplerian velocity distribution is applied to determine the observed line center shift, given by Eq. (E.8). For a massive planet-bearing disk observed under inclination angle i, the Doppler shift of the emission from a single patch of gas in the disk compared to its fundamental frequency ν₀ emerging from a given R,φ point is calculated by $\begin{array}{ccc}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}\mathit{R,\phi ,i}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\nu }}_{\mathrm{0}}}{\mathit{c}}\{{\mathit{u}}_{\mathrm{R}}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\left[\sin \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{+}\cos \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\right]\\ & \mathrm{+}& {\mathit{u}}_{\mathit{\phi }}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\left[\cos \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{-}\sin \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\right]\}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{,}\end{array}$(4)where the radial u_R(R,φ) and azimuthal u_φ(R,φ) velocity components of gas parcels are provided by the hydrodynamic simulations.

The perturbations in the surface density distribution of the disk interior and atmosphere were not taken into account, because the disk atmosphere density given by Eq. (D.2) is independent of the density distribution of the disk even in the gap, as long as there is enough dust in the gap to keep it optically thick. To test the optical thick assumptions, let us for a moment neglect the effect of disk shelf-shadowing that occurres in the inner edge of the gap. Applying the dust density in the disk Σ_d(R = 1   AU) ≃ Σ₀X_d, where X_d = 0.01 is the dust-to-gas ratio and the dust opacity κ_V = 2320   cm²/g, in the gap, depleted to 0.1% of the surrounding density, the optical depth at visual wavelengths along the incident stellar irradiation is τ_V = Σ_d(R = 1   AU)κ_V/δ(R = 1   AU) ≃ 1000. To argue that the shelf shadowing can be neglected we calculated the height of the disk atmosphere and temperature distribution in a planet-free and planet-bearing disk by a 2-D Monte-Carlo code RADMC for dust continuum radiative transfer (Dullemond & Dominik 2004) with the same dust prescription as in our synthetic spectral model. In the planet-bearing disk model the gap is taken to be extended from 0.5 AU to 2 AU and artificially depleted in density to 0.1% that of the planet-free case. The amount of gap depletion is an average value measured in our hydrodynamic calculations. The results are shown in Fig. 2. It is evident that although the height of the disk atmosphere (τ_V = 1 along the line of sight of stellar irradiation) above the mid-plane is decreased in the gap, the temperature in the disk atmosphere is not substantially changed. Because the temperature is not significantly decreased we neglect the surface density perturbations. Note that the effect of the distance of the gap from the central star on the shelf shadowing is not considered. The dust and gas temperature decoupling was also not taken into account. Because the gas has a temperature well in excess of the dust in depleted regions, where the gas column density is  ≪ 10²² (Glassgold et al. 2004; Kamp & Dullemond 2004), we can expect increased contribution to the CO line flux originating from the gap. Contrary to this, if the base of the gap is shadowed completely by the inner wall, the gap does not contribute to the CO emission. Thus a gap could cause a strong permanent line profile asymmetry due to its asymmetric geometry (see Fig. 3), which requires further investigation.

Because the accretion heating has been taken into account according to Eq. (A.9), we had to set an appropriate accretion rate. The accretion rate measured in hydrodynamical simulations gives a value of about 2 × 10^-8   M_⊙/yr. Note that according to Fang et al. (2009) the accretion rate inferred from Hα emission luminosity in the young star population of Lynds 1630N and 1641 clouds in the Orion GMC with an age about 2−3   Myr, is between 10^-10−10^-8   M_⊙/yr. In order to be consistent with the latter, the accretion rate is taken to be 5 × 10^-9   M_⊙/yr for all our models. However, note that at this accretion rate the contribution of the viscous dissipation to the disk continuum is weak in the near-IR band.

4. Results

4.1. Effect of massive planets on the disk structure

Three groups of models were computed through 2000 planetary orbits, with stellar and planetary parameters listed in Table 1. In each run the planet was kept fixed on circular orbit during the first 1000 orbits. After the 1000th orbit the planet was released. It thus felt the (gravitational) backreaction of the disk, resulting in its inward migration. In a good accordance with the expectations, we found that a more massive planet opened a broader and deeper gap. The depletion is 0.25%–0.1% for planets with mass in the range of 1   M_J − 8   M_J. After the first 1000 orbits a quasi steady state disk structure has developed in all simulations, which slightly changed during the following 1000 orbits, when the planet orbit was not fixed anymore.

To shed light on how the giant planet distorts the originally circularly Keplerian gas flow, several snapshots of density and velocity distributions were taken during the 2000th planetary orbit. Figure 3 shows four snapshots of the density and the radial velocity component of the orbital velocity of gas for an 8M_J mass planet orbiting an 1   M_⊙ star (model #8) during one orbit. It is evident that in a disk with an embedded massive planet the overall orbits of gas parcels are non-circularly Keplerian because the radial components of their orbital velocity distribution are strongly departing from zero. Moreover, we found that the elliptic gap, while preserving its shape, precesses slowly retrograde with a period of about 150 planetary orbits.

As we already mentioned, Kley & Dirksen (2006) found that an originally circular disk with an embedded giant planet can reach an eccentric equilibrium state if the mass of the embedded planet is larger than a certain limit (3   M_J for a 1   M_⊙ star). Kley & Dirksen found that for sufficiently wide gaps the growth of the eccentricity is induced by the interaction of the planet’s gravitational potential with the disks material at the radial location of the 1:3 (outer) Lindblad resonance. For smaller planetary masses, this effect is damped mainly by the co-orbital and the 1:2 Lindblad resonances. If the gap is deep and wide enough, which is the case for a giant planet, the above resonances cannot damp the eccentricity-exciting effect appearing at the 1:3 Lindblad resonance, and the disk becomes eccentric. For a more detailed explanation of disk eccentricity growth see Lubow (1991) and Kley & Dirksen (2006). The eccentric state of the disk in our simulations can also be clearly seen in Fig. 3 where the shape of the outer rim of the gap becomes elliptic in the surface density distribution.

The departure of the velocities of the gas parcels from the pure circularly Keplerian circular revolution can be characterized by calculating their eccentricities. Considering the eccentric equilibrium state of the disk, we expect that each gas parcel would move on a non-circular orbit, which may be characterized most conveniently by an average eccentricity value. Therefore, for each disk radius between 0.2 AU and 5 AU, the azimuthally averaged eccentricities of the orbits of the gas parcels were calculated in the following way $\mathit{e}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{h}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\mathit{c}\mathrm{(}\mathit{R,\phi }{\mathrm{)}}^{\mathrm{2}}}\mathrm{d}\mathit{\phi }\mathit{.}$(5)In Eq. (5) c(R,φ) and h(R,φ) stand for $\mathit{c}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\mathrm{=}\mathit{x}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}{\mathit{u}}_{\mathrm{y}}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\mathrm{-}\mathit{y}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}{\mathit{u}}_{\mathrm{x}}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}$(6)and $\mathit{h}\mathrm{\left(}\mathit{R,\phi }\mathrm{\right)}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{x}}\mathrm{(}\mathit{R,\phi }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}_{\mathrm{y}}\mathrm{(}\mathit{R,\phi }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{\mathit{x}\mathrm{(}\mathit{R,\phi }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{y}\mathrm{(}\mathit{R,\phi }{\mathrm{)}}^{\mathrm{2}}}}\mathit{,}$(7)

where u_x(R,φ), u_y(R,φ) and x(R,φ), y(R,φ) are the Cartesian velocity components and coordinates at point R,φ. We found that the averaged eccentricities differ considerably from zero, and each eccentricity curve reaches its maximum near the outer boundary of the gap (Fig. 4). A priori one would expect that the averaged eccentricity values are higher for more massive planets. Indeed we found that the maximum of the azimuthally averaged eccentricity is monotonically increasing with a planetary mass in the range of 1   M_J ≤ m_pl ≤ 8   M_J. The eccentricity curve for the 8   M_J mass planet peaks about e_max ~ 0.3, while the peak stays well bellow 0.1 for a 1   M_J mass planet. Note that a very similar behavior of the disk eccentricity has already been found by Kley & Dirksen (2006). In their cases, however, the maximum of the eccentricity curves is somewhat lower than in our cases, and at least 3   M_J is required to set the disk into eccentric state. This can be the consequence because contrary to Kley & Dirksen (2006), we used an α-type viscosity in our simulations. The higher eccentricities we found are plausible inasmuch as Kley & Dirksen (2006) found that the eccentricity of the disk is increasing with decreasing viscosity and in our approach the kinematic viscosity ν(R) = αH²Ω_K(R) measured in dimensionless units is 2.5 × 10^-6 at R = 1, which is smaller than the ν = 1 × 10^-5 used by Kley & Dirksen (2006).

As one can see, a giant planet has substantial impact on the density and velocity distributions of its host disk. It is an essential question, whether the velocity perturbations appear in the line-of-sight velocity with substantial strength. Figure 5 shows the subtracted distribution of the line-of-sight velocity in the perturbed and circularly Keplerian case. It is evident that the line-of-sight velocities show a significant departure from the circularly Keplerian fashion because the difference is non-vanishing. The radial velocity component of the orbital velocity and more importantly the line-of-sight velocity distributions have variable patterns following the planet on the top of a permanent excess seen at 0° and 180° position angle. Note that the the permanent pattern precesses slowly (with  ~150 orbital periods), retrograde to the planet. The disk inclination angle i is taken to be 40° in the calculation of Fig. 5 and the disk was rotated to the line-of-sight in a way that the line-of-sight velocity component of dynamically perturbed gas parcels is maximized. Consequently we had expected that not only significant distortions appear in the line profiles, but that they vary in time within the orbital time scale of the giant planet. Because the deviation from circularly Keplerian velocity is in the range of −7.92−4.88   kms^-1, the width of variable component in the line profile should be  ~10   kms^-1, depending on the inclination angle.

4.2. Distortion of CO lines

Below we address the following questions: (i) Does the line-of-sight component of the non-circularly Keplerian velocity distribution (Fig. 5) result in significant distortions in the CO spectral line profiles? (ii) Are these distortions varying in the planets orbital timescale? (iii) How do the distortions depend on the inclination angle, the planetary and stellar mass, the orbital distance of the planet, the size of the inner cavity and the disk geometry (if the disk is flared or not)? (iv) Are these distortions strong enough to be detected by high-resolution near-IR spectral measurements? A complete set of answers to the above questions may give us a novel method to detect massive planets embedded in protoplanetary accretion disks.

As was shown in Fig. 4 the disk eccentricity is considerable near the gap in all planet-bearing disk models. If the gas parcels move on pure elliptic orbits in the gap, the CO line profile becomes permanently asymmetric. This effect is illustrated in Fig. 6, where we have calculated the emerging line profile of a disk with an eccentricity e = 0.2 in the gap between 0.8   AU−2.8   AU hosted by a 1   M_⊙ star viewed 20°, 40° and 60° degree of inclination angles. Thus we can expect that the non-circularly Keplerian property of the velocity field shown in Figs. 3 and 5 will break the symmetry of the spectral lines. Although the orbits of the gas parcels can be characterized by average eccentricities, their orbits are more complicated than those of a pure elliptic motion, because of the perturbations induced by the planet.

First we examined the CO emission line profile distortions from a disk in which a 8 M_J planet orbits a 1 M_⊙ star (model #8). As expected, the strongly non-circularly Keplerian velocity flow significantly modified the CO line profiles. In Fig. 8a we display the V = 1–0 P(10) line profiles for three different inclination angles. The line profiles show a strongly asymmetric double-peaked shape. Studying Fig. 8a two distortion patterns of the spectral line profiles can be recognized. An excess can be seen near the red peak of the lines at about 11   kms^-1, 21kms^-1 and 29   kms^-1, for inclinations i = 20°, 40° and 60°, respectively. Moreover, a deficiency appears around the blue peak of each spectral line at slightly lower velocities. The above velocities correspond to the orbital velocity of the gas revolving around a 1   M_⊙ star at 0.8–0.9 AU, i.e. near the inner boundary of the gap, taking into account the appropriate inclination angles. Smaller-scale distortions also appear in the red and blue peaks as well, corresponding to the 0.8–2.9 AU region, i.e. the whole gap.

Calculating the line profiles at different orbital phases during one orbit, we found that the shape of the already distorted line profile varies in time, yielding a clear dependence on the orbital position of the planet, see results in Fig. 9 for model #8. The permanent asymmetry in lines is caused by the permanent velocity pattern seen at PA 0° and 180° in Fig. 5. Note that the expression “permanent” is not entirely correct because the velocity pattern slowly (~150 orbital periods of planet) precesses retrograde to the planet, but we still use it henceforward. The variable component indicated with arrows is moving in regions between  ~ ±25   kms^-1 due to the variable pattern following the planet. Note that the maximum width of variable component is  ~10   kms^-1, as is expected, which is resolvable by CRIRES, whose maximum resolution is 3   kms^-1. In this particular model the time scale of variations is on the order of weeks, because the orbital period of the planet is one year, thus the time elapsed between snapshots is approximately 18 days.

4.3. Planetary and stellar masses

To gain further insights, it is useful to study the influence of the masses on the strength of the spectral line distortions. First we have investigated the cases in which a 1   M_⊙ star hosts a 1   M_J, 3   M_J, 5   M_J, and 8   M_J mass planet at 1 AU, corresponding to models #5-#8. Parts of the results are presented in Fig. 8b for model #5 (m_pl = 1   M_J, m_∗ = 1   M_⊙) and in Figs. 8a and 9 for model #8. One can conclude that for a less massive planetary companion the line profile distortions are weaker. The same conclusion can be obtained by the growing influence of the increasing planetary mass on the disk eccentricity, see Fig. 4.

A natural way to study the effect of the mass of the central star on the line profile distortions would be that the planetary mass is kept fixed, and the stellar mass is increased. On the other hand, we recall that dimensionless units were used in hydrodynamical simulations thus the mass of the planet is expressed in stellar mass units. Thus changing either the planet mass or the stellar mass is equivalent to changing the ratio between the planetary and the stellar mass. Therefore, one could conclude from the previous results that the larger the stellar mass, the stronger the line profile distortions. The situation is however a bit more complicated; if we change, for instance, the stellar mass, the luminosity of the star changes, which has an influence on the temperature distribution in the disk atmosphere and interior as well. The effect of the stellar mass can be therefore studied if the planetary-to-stellar mass ratio is kept fixed, while the mass along with the surface temperature and the stellar radius is changed in the synthetic spectral model. Due to the effect of the increased/decreased flux of stellar irradiation in case of the larger/smaller stellar mass, the outer boundary of the CO emitting region will be moved farther/closer to the star, and increased/decreased in size. To investigate the dependence of the observable line profile distortions on the stellar mass, we have recalculated the above presented line profiles for a 0.5   M_⊙ (models #1–#4) and a 1.5   M_⊙ (models #9-#12) central stars. Part of the results are shown in Fig. 8c (for model #4) and Fig. 8d (for model #12), respectively. Comparing the line profiles to those presented in Fig. 8a (for model #8), it is evident that the profiles are considerable narrower for 0.5   M_⊙ and broader for 1.5   M_⊙ models, due to the change in the Keplerian angular velocity (Ω_K(R) = (Gm_∗/R³)^1/2). Moreover, the line-to-continuum ratio at the maximum is slightly decreased and increased about a same amount for disks with 0.5 and 1.5   M_⊙ star. The change in line-to-continuum ratio can be easily explained: for a smaller mass star, the disk temperature is also decreased due to the decrease in stellar luminosity in a way that the ratio of the CO line to the disk plus star continuum is decreased too. The opposite is true for larger stellar masses. Finally, we can conclude that for a given planet-to-star mass ratio the larger the central stellar mass, the larger the observable CO line profile distortion caused by the giant planet. Here we have to note that although the observational data on T Tauri stars (Akeson et al. 2005) show that the size of the disk inner cavity is increasing with increasing stellar luminosity (which can be the consequence of the increased dust sublimation radius, see Eq. (3)), resulting in smaller continuum and stronger relative line strength, in our cases the size of inner cavity stayed fixed in order to make models comparable.

Analyzing Fig. 8b we conclude that the line profile distortions in model #5 (m_pl = 1   M_J, m_∗ = 1   M_⊙, 1 AU orbital distance) are so small that the giant planet signal below 1   M_J is strongly suppressed. Note that a disk with a 1   M_J mass planet embedded into it will not get into the eccentricity state at all, see Fig. 4. Planets orbiting the host star closer than 1 AU, however, can induce stronger distortions, as shown in the following section.

4.4. Orbital distance of the planet

In this section we investigate the effect of the orbital distance of the planet to the strength of the spectral line distortion. One would initially expect that the spectral line distortions are stronger/weaker as the orbital distance of the planet decreases/increases. We recalculated the hydrodynamical simulations in “tight” and “wide” models, where the planets orbit at 0.5 and 2 AU distances from the central star, respectively. In good agreement with the expectations, we found that the planet signature in the spectral line was strongly suppressed in the “wide” versions of models #1–#8, where m_∗ = 0.5   M_⊙ and m_∗ = 1   M_⊙. The reason for this is obvious: in wide systems the planet is orbiting in cold regions (the atmosphere temperature is varied between 200–350 K at 2 AU, depending on the stellar mass), where the atmospheric CO emission measured to the continuum is weak. As can be seen in Fig. 8e, the “wide” version of model #8 (m_pl = 8   M_J, m_∗ = 1   M_⊙, 2 AU orbital distance) shows considerably weaker planet signatures in the CO line profile than the original 1 AU model, see Fig. 8a for comparison. The distortion patterns appear at lower velocities owing to lower orbital velocities in 0.5   M_⊙ models. On the other hand, we found about the same strength of the planet signal in the “wide” version of model #12 (m_pl = 12   M_J, m_∗ = 1.5   M_⊙, 2 AU orbital distance) as in model #8. This is expected, since in the latter case the luminosity of the host star is high enough to produce sufficient CO excitation even at 2 AU.

The resulting line profile calculated in the “tight” version of the model #8 (m_pl = 8   M_J, m_∗ = 1   M_⊙, 0.5 AU orbital distance) is shown in Fig. 8f. Because the giant planet orbits closer to the host star, where the disk atmosphere is heated to higher temperatures (the atmosphere temperature varies between 450–700 K at 0.5 AU, depending on the stellar mass), the emission emerging from the perturbed regions is stronger than in the original 1 AU models. It is also evident that the distorted peaks are slightly shifted to higher velocities compared to the original 1 AU models, because the orbital velocity of the gas parcels in dynamically perturbed orbits is higher in “tight” models. The disk is more eccentric on average in the latter case (Fig. 7), but the strength of the line profile distortion is increasing as the planetary orbital distance decreases.

In order to demonstrate that smaller mass planets can also be detected in closer orbits we have calculated the line profile in a “tight” version of model #2, m_pl = 1.5   M_J, m_∗ = 0.5   M_⊙, 0.5 AU orbital distance. We find that one could observe roughly the same strength line profile distortions as in model #8 (m_pl = 8   M_J, m_∗ = 1   M_⊙, 1 AU orbital distance).

4.5. Inner cavity size

Now we turn our attention to the question under which circumstances planets can be detected orbiting farther away from the star (farther than 2 AU for instance). In Sect. 3.2 we argued that there is a region around the star called inner cavity from which the CO is already depleted. In our simulations presented so far, we assumed that the radius of the inner cavity is fixed at 0.2 AU irrespective of the stellar luminosity. On the other hand, the size of the inner cavity is presumably dependent on the stellar mass, and indeed this is confirmed by observations (Akeson et al. 2005). As we mentioned in Sect. 3.2, the size of the inner cavity has a significant influence on the overall line-to-continuum ratio. In disks dynamically perturbed by planets orbiting at 1 AU and with a significant amount of CO lying inside 0.2 AU, the majority of CO flux is emerging from the unperturbed innermost regions, which eventually could smear out the planet signatures. On the contrary, if the size of the inner cavity is larger, the contribution of the perturbed regions to the total CO flux becomes stronger. In Fig. 8g we present the line profiles obtained in model #8, in which the inner cavity is increased to 0.4 AU in size, while the 8   M_J mass planet was still orbiting at 1 AU. As can be seen, the overall line-to-continuum ratio is decreased due to the absence of hot gas, but the level of the asymmetric pattern profile is more apparent than in models with an inner cavity extending only to 0.2 AU, see Fig. 8a for comparison. For larger inner cavities the continuum level of the disk interior and inner rim are decreased too, resulting in a moderate weakening of the line-to-continuum ratio and slight strengthening of planet signal. In this case the detection possibility of a planet orbiting at larger distances (> 2   AU) is growing.

4.6. Disk geometry

It is well known that many T Tauri sources present flatter than λF_λ ~ λ ^−4/3 SEDs. One attempt to explain this SED flattening was to assume that the disk is flaring (Kenyon & Hartmann 1987), resulting in a geometrically thicker disk atmosphere H(R) ~ R^γ + 1, where γ is the flaring index. In the approach of Chiang & Goldreich (1997) the flaring index is γ = 2/7. In order to investigate the effect of disk geometry on the strength of the line profile distortions we recalculated several models using flat thermal disk approximation (the hydrodynamical inputs were the same). An example of the flat version of model #8 (m_pl = 8   M_J, m_∗ = 1   M_⊙, 1 AU orbital distance) is shown in Fig. 8h. Note that in order to have the same peak line-to-continuum ratio in flat as in flared models, the amount of emitting CO is arbitrary increased by decreasing the dust-to-gas ratio to 1.2 × 10^-3 and keeping the disk mass unchanged. The reason for this is that although the atmospheric temperature does not differ in flared and flat disk models (see the atmospheric temperature given by Eq. (A.5)), the atmospheric surface density of CO is substantially smaller in flat than flared disk models (see the dependence of surface density given by Eq. (D.2) on the grazing angle, given by the first term of Eq. (A.2) in flat disk models). It is evident that the total line-to-continuum ratio is higher in the flat than in the flared disk model, but the line width at the foot of the line profile does not change. According to our calculations the line profile distortions in flat models compared to the flared one are less significant in “wide” systems where planets orbit at larger distances. This can be explained by the effect of flaring getting stronger with increasing distance to the host star.

4.7. Observational considerations

To observe the planet-induced distortions of the CO line profiles that we simulated, a spectroscopic facility is required that meets two basic demands: it needs sufficient spectral resolution to properly resolve the sub-structure in the line profile and sufficient sensitivity to bring out the low-contrast line profile distortions. In this section, we investigate the observability of the modeled effect in the context of contemporary facilities, and give a brief outlook into the ELT era. We present a quantitative example for the Cryogenic High Resolution Echelle Spectrograph (CRIRES) at the VLT (Kaeufl et al. 2004), which is currently the most powerful high-resolution spectrograph covering the M-band around 4.8   μm, capable of observing CO ro-vibrational spectra.

The main limiting factor for ground-based thermal infrared observations is the Earth atmosphere, which strongly absorbs the radiation from astronomical sources and causes a high thermal background. Particularly strong telluric absorption and emission is present at the wavelengths of low-excitation CO transitions. Within  ~10   kms^-1 from the line center the telluric transmission is so low and the background emission so high that these spectral regions are not accessible. However, depending on the location of the source with respect to the ecliptic, the CO features are Doppler shifted by up to  ± 30   kms^-1 due to the Earth’s orbital motion, making the whole CO line accessible over the course of several months, in principle. Note though that we can never measure the entire CO line profile simultaneously, there will always be an approximately 20   kms^-1 wide “gap” in our spectral coverage at any observing epoch.

In order to test whether we would be able to detect the planet-induced permanent profile asymmetries with currently available instrumentation, we have simulated a VLT/CRIRES observation of model #8. We have assumed the source to have a brightness of M = 6   mag, and calculated the achieved SNR as a function of wavelength using the VLT/CRIRES exposure-time calculator⁴. This calculation takes into account the telluric absorption and emission, the quality of the adaptive optics correction, the system throughput and detector characteristics. We note that our simulation only considers the signal-to-noise ratio, and assumes that systematic effects due to, e.g., time-variable telluric absorption, can be well calibrated. The achieved SNR can be scaled according to $\mathit{SNR}\mathrm{\propto }{\mathrm{10}}^{\mathrm{0.4}\mathit{M}}\sqrt{{\mathit{T}}_{\mathrm{int}}}\mathit{,}$(8)where M denotes the M-band magnitude and T_int the integration time. In Fig. 10 we show the emerging V = 1–0P(10) line profiles for model #8 (upper figure) and the atmospheric transmission (lower figure) for an assumed integration time of 1 h on a source of M = 6   mag. The model spectrum has been convolved with the 3 km s^-1 CRIRES instrumental resolution. We have assumed zero relative velocity between the source and the Earth at the time of the observation. The effects of the telluric CO line are clearly seen around the line center, where no meaningful measurement can be obtained. The line wings, where the effects of the planet are most prominent, can be well measured.

In order to quantify the significance with which a planetary signal (the most prominent permanent asymmetric pattern) may be detected, we have tried to fit the modeled observation with both planet-free and planet-bearing models. In Fig. 11 we show the goodness-of-fits (reduced χ²) as a function of integration time for both cases (with dashed line for planet-free and solid line for planet-bearing disk models) for 20°, 40°, and 60° inclination angles. Above a certain integration time (presented by vertical lines) the planet-bearing models yield a more than 50% better reduced χ² fit to the simulated data than those of the planet-free models. The asymmetric pattern induced by an embedded planet increase and thus become easier to detect with disk inclination.

Because the fundamental band ro-vibration R(0), R(1), R(12), R(13), R(19)-(21), R(25), R(26) and R(30)-(33), P(6-12) and P(22-32) CO lines are not blended by the higher vibrational lines⁵, the line profiles can be averaged, using appropriate scaling, resulting in an increased signal-to-noise ratio. Note that the unblended lines with higher J > 25 or lower J < 2 rotational quantum numbers are too weak for averaging: using them would introduce only additional noise to the averaged profiles. Also note that, since the higher vibrational levels are not excited at all in disk atmospheres similar to that used in our models, the blended lines can also be used for averaging, resulting in an even more increased signal-to-noise ratio.

We conclude that in the confines of our model, the permanent asymmetric signals produced by an 8   M_J giant planet orbiting at  ~1   AU embedded in the disk of a T Tauri star of brightness M = 6   mag would be detectable with VLT/CRIRES with an integration time of 1 h. Moreover, because the variable component at the top of the permanent asymmetry has a width of  ~5−10   kms^-1, exceeding the CRIRES  ~3   kms^-1 resolution, and  ~10% magnitude, there is a certain possibility to strengthen its planetary origin.

With the next generation of giant, ground-base telescopes our sensitivity to low-contrast features in infrared spectra will increase greatly. We have used the ESO exposure-time calculator for the ELT (version 2.14) to simulate an observation with the proposed first generation infrared instrument Mid-infrared E-ELT Imager and Spectrograph (METIS) (Brandl et al. 2008). The exposure-time calculator necessarily incorporates planned telescope and instrument specifications, because this facility is yet to be built. We find that ELT/METIS yields an increase in sensitivity of a factor of  ~30 compared to VLT/CRIRES. Thus also sources that are substantially fainter or have weaker spectral features than those modeled here in the context of VLT/CRIRES are within reach of the observational facilities of the next decade.

5. Conclusions

We have shown that our semi-analytical 2D double-layer thermal disk model (in which the disk layers are heated mainly by stellar irradiation, and the emission of the disk inner rim is accounted for) is capable to reproduce the double-peaked Keplerian line profiles in CO fundamental ro-vibrational band in young (about 2.5   Myr old) T Tauri-type protoplanetary disks assuming canonical dust and gas properties. However, note that several T Tauri stars show centrally peaked symmetric profiles instead, which can be attributed to an enhanced CO flux formed in distant regions (up to 10 AU), presumably above the disk atmosphere. This can be explained by non-LTE heating, such as UV fluorescence (Krotkov et al. 1980), or dust gas temperature decoupling (Glassgold et al. 2004; Kamp & Dullemond 2004), resulting in highly excited CO emissions. These effects were not taken into account in our calculations. Thus the gas temperature at disk atmosphere provided by our thermal model is presumably somewhat lower than the realistic one. Our predictions, however, demonstrate a “the worst case” because the planet signal is getting stronger with increasing atmospheric temperature. It is revealed that significant line profile distortions appear in the CO fundamental ro-vibrational band caused by a giant planet. As we have demonstrated, a giant planet with a mass of at least 1   M_J can perturb the disk sufficiently to be detected by permanent asymmetry in line profiles and its short timescale variations. The permanent asymmetry in line profiles can be interpreted by the disk being in an eccentric state only in the gap, but the short timescale variability is certainly connected to the local dynamical perturbations of the orbiting planet. Depending on the signal-to-noise ratio and the resolution of the spectra, even lower-mass giant planets can be detected revolving on close orbits in disks observed with high-inclination angles and with large size inner cavities. The summary of our findings is:

• 1.

  The dynamical perturbation induced by a giant planetsignificantly distorts the CO ro-vibrational fundamentalemission line profiles.

• 2.

  The line profile distortions are more apparent for larger inclination angles.

• 3.

  The main component of line profile distortions is a permanent asymmetry. Its position in the line profile depends on the stellar mass and orbital distance of the planet. The permanent asymmetry shape depends on the orientation of the gas elliptic orbit with respect to the line of sight.

• 4.

  The line profiles are changing with time, correlated to the orbital phase of the giant planet. The timescale of the variation is on the order of weeks, depending on the orbital period of the planet, i.e, the orbital distance and the mass of the host star. The magnitude of variable component is  ~10%, its width is  ~10   kms^-1, depending on the inclination angle.

• 5.

  The planet signal becomes stronger with increasing planetary mass.

• 6.

  The planet signal strengthens/weakens with increasing/decreasing mass of the host star, neglecting the effect of the stellar luminosity on the size of the inner cavity.

• 7.

  The planet signal strengthens/weakens with decreasing/increasing orbital distance of the planet.

• 8.

  The size of the inner cavity has a strong influence on the giant planet observability. If the size of the inner cavity is substantially smaller than 0.2 AU, we can detect a giant planet only at close (about 0.5 AU) orbit. On the contrary, if the cavity is larger (0.4 AU), a distant planet (in 2 AU orbital distance) can also be detected.

• 9.

  The influence of disk geometry on the planet signal is modest. With a decreasing flaring index the strength of planet signal does not change although the line profiles are strengthened, except in “wide” flat models where the planet signal substantially suppressed.

• 10.

  The lowest mass of giant planet orbiting at 1 AU that still can be detected is 0.5   M_J for a 0.5   M_⊙ star in the confines of our models.

• 11.

  By high-resolution near-IR spectroscopic monitoring with VLT/CRIRES, giant planets with  ≥ 1   M_J mass orbiting within 0.2–3 AU in a young disks may be detectable, if other phenomena do not confuse, mimic, or obscure the planet signature. The level of confidence is growing with increasing inclination angle.

In the light of our findings, we propose to observe T Tauri stars that shows an appreciable amount of CO in emission with high-resolution near-IR spectrograph to search for giant planet companions. The asymmetric and time-varying double-peaked line profiles can be explained by strongly non-circularly Keplerian gas flow in the disk caused by the giant planets, as was presented. Although higher V ≥ 2 vibrational states of CO are not excited in our models, overlaps can occur, resulting in asymmetric line profile. To avoid this, one should consider only non-blended lines for profile averaging such as P(6)-P(12), and R(0), R(1), R(12), R(13), R(19)-(21), R(25), R(26) and R(30)-(33). Another possible source of permanent line profile distortions could be a large companion deforming the circumprimary disk to fully eccentric in mid-separation binaries (e.g. Kley et al. (2008); Regaly et al. in prep.). Last but not least, the distorted PSF, caused by tracking errors of the telescope or unstable active optics during an exposure, can induce artificial signals. A real signal from a bipolar structure will change sign when observed using antiparallel slit position angles whereas the artificial signal will not. Thus artificial signatures can be successfully identified with the comparison of the spectra of disks obtained at two antiparallel slit position angles, e.g. 0° and 180° (Brannigan et al. 2006). By detecting short time-scale (weeks to months) variations of the asymmetric profiles with similar patterns as presented in Fig. 9, we could infer the presence of a giant planet in formation.

The proposed instrument METIS for the ELT (Brandl et al. 2008), providing high-resolution (λ/Δλ = 10⁵) infrared spectroscopy with an increase in sensitivity of a factor of  ~30 to that of VLT/CRIRES will give the opportunity to even decrease more the mass and brightness limits of possible giant planets. If we revealed for instance the distance of the observed planets to the star as a function of the age of the star-disk system, we would be able to constrain planetary migration models. Moreover, determining the youngest star host planet-bearing disks could also provide constraints on planet-formation models.

Although our assumptions are reasonable from the point of view of finding a simple model of the signal of dynamical perturbations, it can be developed further of course. For example, if the gap is cleared by the planet to a fraction of the surrounding density, the disk is not directly irradiated, but is shadowed by the inner wall of the gap (Varnière et al. 2006). This mild heating causes a temperature drop at the inner gap wall, which results in a depressed CO flux. On the other hand, the outer edge of the gap is directly illuminated by the star, resulting in a slight temperature increase, which produces stronger CO emission from the regions where the eccentricity reaches its maximum (see e.g. Fig. 4). According to our results the gap is strongly non-axisymmetric, thus 3D radiative transfer using 3D density structure needed to calculate the real temperature profile in perturbed disks. Moreover, as the gap is depleted in large grains because of the trapping of dust grains larger than  ~10   μm (Rice et al. 2006), the CO flux normalized to the continuum is strengthened. The gas and dust can be thermally uncoupled (see e.g., Glassgold et al. (2004); Kamp & Dullemond (2004); Woitke et al. (2009)) in the tenuous gap, which may result in CO being overheated to the dust. These effects could cause additional permanent line profile distortions due to the elliptic shape of the gap seen in density distribution (Fig. 3). If we consider an excess emission originating from the planetary accretion flow suggested by Clarke & Armitage (2003), or the density waves in disk surface caused by the orbiting planet, it would cause a varying planet signal. To these caveats, one might also consider the disk turbulence (e.g., due to MRI effects), which may add random asymmetries and “blobbiness” to the line profiles. We will investigate these effects in a forthcoming paper using a more sophisticated thermal disk model.

CO exists inside the dust sublimation radius in some T Tauri stars (Carr 2007) or Herbig Ae/Be stars (Brittain et al. 2009), where the disk is optically thin in the absence of dust. One possibility to excite CO in the optically thin inner disk close to the stellar surface is the UV fluorescence, which is far from LTE (Krotkov et al. 1980) applied in our calculations. Depending on the efficiency of X-ray heating the gas temperature can be as high as  ~4000−5000   K. Somewhat below this temperature, where CO still exists, the high-excitation vibrational levels of CO such as V ≥ 2 are populated, resulting in significant V = 2 → 1, V = 3 → 2, etc. ro-vibrational lines. According to Skrutskie et al. (1990), the inner cavity in transitional disk may be formed by a close giant planet. Indeed, Lubow & D’Angelo (2006) found that an 1−5   M_J mass planet significantly lowers the accretion of dust and gas to 10%–25% of the accretion rate outside the gap. We expect that the planetary signal forming in the optically thin inner disk is also detectable, which is strongly supported by the revealed line profile distortion strengthening with decreasing orbital distance of the planet. The UV pumping requires considerable UV flux from the central star, which is a characteristic of Herbig Ae/Be stars rather than smaller T Tauri stars. Considering that the CO may be excited up to to 10 AU by UV fluorescence (Brittain et al. 2007, 2009) in Hebig Ae/Be stars, planet orbiting at larger distances than 2 AU may be also detectable by CO ro-vibrational line profile distortions.

Spectroscopy of T Tauri stars shows emission of molecules such as H₂O, OH, HCN, C₂H₂, and CO₂, as well. Nevertheless CO is more abundant than these molecules by  ~10 times, only H₂O could reach the abundance of CO, predicted by recent models that calculate the vertical chemical structure of the gas in disk atmosphere (e.g. Glassgold et al. 2004; Kamp & Dullemond 2004; Woitke et al. (2009)). Indeed in some cases, like AA Tau (Carr & Najita 2008), and AS 205A and DR Tau (Salyk et al. 2008), the rotational transitions of H₂O dominate the mid-infrared (10−20   μm) spectra, suggesting that H₂O is abundant in disk atmospheres. The strong water emission could be the consequence of turbulent mixing that carries molecules from the disk midplane, where they are abundant, to the disk atmosphere (Carr & Najita 2008), or the effects of an enhanced mechanical heating of the atmosphere (Glassgold et al. 2009). While the H₂O ro-vibrational lines probe the inner regions of disk out to radii  ≥ 2   AU, the rotational lines are produced between radii 10-100 AU (Meijerink et al. 2008). The rotational and ro-vibrational line profiles of water are also subject to distortions caused by the dynamical perturbations of a giant planet. Thus it is reasonable to search for signature of giant planets in the high-resolution spectra of H₂O in the mid-infrared band (for planets orbiting int the inner disk), and rotational spectra of H₂O in the far-infrared band (for planets orbiting in the outer disk), as well. Because H₂O is heated by stellar X-rays and sub-thermally populated beyond  0.3 AU, X-ray-heating and non-LTE level population treatment is needed to calculate water lines (e.g., Meijerink et al. 2008; Kamp et al. 2010).

Finally, we point out some noticeable resemblance to our findings revealed in observed CO line profiles. The V836 Tau transitional disk shows strongly distorted 4.7   μm CO ro-vibrational line profile presented by Najita et al. (2008). The averaged CO line profiles shows very similar features to our results. According to Najita et al. (2008) the possibility of a massive planet is also restricted by current limits on the stellar radial velocity of V836 Tau, which constrain the mass of a companion within 0.4−1   AU to 5−10   M_J. Because the disk extents to a limited range of radii (~0.05−0.4   AU), the dynamical perturbation of a close orbiting planet with a mass smaller than 5   M_J presumably could cause the observed line profile distortions. Note that radial velocity measurements indicate that no planet is present larger than 1−2   M_J closer than  ~0.5   AU (Prato et al. 2008). The transitional disks that encompass the young sources SR 21 and HD 135344 B have been observed by VLT/CRIRES (Pontoppidan et al. 2008) and they already show clear line profile asymmetries in CO fundamental ro-vibrational band. Our scenario of a planet is a reasonable explanation of the asymmetry. The presence of a planet was also proposed by Grady et al. (2009) for HD 135344 B and Eisner et al. (2009) for SR 21, based on SED and visibility data, respectively. EXLupi (prototype of EXor type young variable stars) showed remarkable line profile variations during its 2008 outburst (Goto et al. 2010, in prep.). Goto et al. revealed that the CO emission has a spatially multiple origin. The quiescent component forms in the outer optically thick disk, while the outburst component in the inner optically thin gas disk, where the higher vibrational levels of CO are presumably excited by UV pumping. These high-excitation V = 2 → 1 and V = 3 → 2 lines are double peaked, strongly asymmetric and variable on a short (weekly) time scale. Thus for example, by measuring short timescale (on the order of weeks or less) periodic variations in the CO ro-vibrational spectra of these sources, by mid-IR monitoring observations would eventually indirectly detect an embedded Jupiter-like planet in birth.

---

Acknowledgments

It is a pleasure to acknowledge helpful discussions with Ewine van Dischoeck, Klaus Pontoppidan and Attila Juhász. This research has bee supported in part by DAAD-PPP mobility grant P-MÖB/841/ and “Lendület” Young Researcher Program of the HAS. We are grateful to L. Kiss who helped us to substantially clarify the text. We also thank the anonymous referee for thoughtful comments that helped to significantly improve the quality of the paper.

References

Appendix A: Temperature distribution in the disk

For the sake of completeness we present here the derivation of the dust temperature distributions in the disk interior and atmosphere according to Chiang & Goldreich (1997). Below we use the local thermodynamical equilibrium (LTE) assumption everywhere. The stellar flux in the disk atmosphere at a distance R from the stellar surface is ${\mathit{F}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\simeq }\frac{\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathrm{2}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}\mathit{,}$(A.1)where δ(R) is the grazing angle of the incident irradiation, R_∗ and T_∗ are the stellar radius and surface temperature, respectively, and σ is the Stefan–Boltzmann constant. Here we assumed that only half of the stellar surface is visible from a given point at a distance R from the star in the disk atmosphere. According to Chiang & Goldreich (1997) the grazing angle can be given by $\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\mathrm{5}}\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)\mathrm{+}\frac{\mathrm{8}}{\mathrm{7}}{\left(\frac{{\mathit{T}}_{\mathrm{\ast }}}{{\mathit{T}}_{\mathrm{g}}}\right)}^{\mathrm{4}\mathit{/}\mathrm{7}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{7}}\mathit{,}$(A.2)in a flared disk assumed to be in hydrostatic equilibrium. In Eq. (A.2) T_g is the gravitational temperature at which the thermal energy of gas parcel balances the gravitational energy at the stellar surface, which can be given by ${\mathit{T}}_{\mathrm{g}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathit{m}}_{\mathrm{H}}}{\mathit{k}{\mathit{R}}_{\mathrm{\ast }}}\mathit{,}$(A.3)where m_H is the mean molecular weight of the gas, G and k is the gravitational constant and Boltzmann constant, respectively. Note that for a non-flared, i.e. flat disk the grazing angle is defined by only the first term of Eq. (A.2).

By definition the optical depth of the disk atmosphere along the incident stellar irradiation is τ_V = 1 in the optical wavelengths. The optical depth of the atmosphere at near-IR heated to T_atm,irr(R) perpendicular to the disk plane is τ_Vδ(R)ϵ_atm(R), where ϵ_atm(R) = (T_atm,irr(R)/T_∗)^β is the dust emissivity, and β is the power law index of the absorption coefficient of the dust, see Appendix C for details. The total flux emitted inward and upward by the optical thin atmosphere can be given by ${\mathit{F}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathit{ϵ}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{\sigma }{\mathit{T}}_{\mathrm{atm}\mathit{,}\mathrm{irr}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathit{.}$(A.4)In LTE the absorbed flux of atmosphere equals to the emitted one (F_∗(R) = F_atm,irr(R)), thus the disk atmosphere is heated to the temperature ${\mathit{T}}_{\mathrm{atm}\mathit{,}\mathrm{irr}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\left(\frac{\mathrm{1}}{\mathrm{4}}\right)}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\beta }\mathrm{)}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{2}\mathit{/}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\beta }\mathrm{)}}{\mathit{T}}_{\mathrm{\ast }}$(A.5)by stellar irradiation.

Assuming that the disk interior with a temperature T_int,irr(R) is optically thick to its radiation everywhere in our computational domain,⁶ the flux F_int,irr(R) emitted by the disk interior is ${\mathit{F}}_{\mathrm{int}\mathit{,}\mathrm{irr}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{irr}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathit{.}$(A.6)Considering that the absorbed stellar flux will be re-emitted by the disk atmosphere upward and downward only half of this flux heats the interior, i.e. 0.5F_atm(R) = 0.5F_∗(R) = F_int,irr(R). Accordingly, the temperature of the disk interior set by the stellar irradiation can be given by ${\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{irr}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\left(\frac{\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathrm{4}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}_{\mathrm{\ast }}\mathrm{·}$(A.7)Owing to the accretion process a significant amount of gravitational potential energy has to be dissipated by a viscous processes. According to Lynden-Bell & Pringle (1974), in a steady state disk⁷ with a constant accretion rate Ṁ, the flux F_acc(R) released by the disk mid-plane due to the change of potential energy can be given by ${\mathit{F}}_{\mathrm{acc}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{Ṁ}}{\mathrm{8}\mathit{\pi }}\left(\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right){\mathit{R}}^{-3}\mathit{.}$(A.8)In this way the disk mid-plane is heated to the temperature ${\mathit{T}}_{\mathrm{acc}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\left[\frac{\mathrm{3}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{Ṁ}}{\mathrm{8}\mathit{\pi \sigma }}\left(\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{R}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{R}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{4}}\mathit{.}$(A.9)Here we have to note that only half of the total accretion power is involved in Eq. (A.8), the remaining is stored in the kinetic energy of orbiting gas. This energy should be radiated away by the disk boundary layer (Popham et al. 1993, 1995) or by the accreting material in the funnel flow formed along the magnetic field lines in magnetospheric accretion model, see Hartmann et al. (1994); Bouvier et al. (2007). Because the disk boundary layer and the funnel flows are confined into such small volumes that the gas temperature reaches about 10 000 K emitting in the UV band, its radiation is not taken into account.

Now let us take into account the heating due to viscous dissipation in disk interior using the superposition principle, i.e. the radiation at a given location is the sum of the radiations corresponding to different heating sources. To determine the disk interior temperature T_irr,acc(R) caused by the viscous dissipation first assume that the flux emitted by the optical thick disk interior is ${\mathit{F}}_{\mathrm{int}\mathit{,}\mathrm{acc}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{acc}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathit{.}$(A.10)Because the disk mid-plane radiates the accretion flux into two directions 0.5F_acc(R) = F_int,acc(R), the disk interior is heated to ${\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{acc}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{T}}_{\mathrm{acc}}\mathrm{(}\mathit{R}\mathrm{)}$(A.11)by accretion. Taking into account the heating of the disk interior by irradiation of atmosphere and viscous dissipation together (simply summing the radiation fluxes), the resulting temperature of disk interior is ${\mathit{T}}_{\mathrm{int}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\left(}{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{irr}}\mathrm{\right(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{acc}}\mathrm{\left(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{\right)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}$(A.12)

Appendix B: Emission from the disk inner rim

To determine the temperature profile of the disk interior close to the disk inner edge, where the simple double layer assumption cannot be applied, we first assume that the rim is a perfect vertical wall. Any given optically thick vertical slab with a thickness of dR at R + dR distance from the central star is irradiated by the neighboring hotter slab located at R. Assuming that half of its radiation is received by the one at R + dR, the emitted and irradiated fluxes are in equilibrium in LTE, i.e. $\frac{\mathrm{1}}{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{rim}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{=}\mathit{\sigma ϵ}\mathrm{(}\mathit{R}\mathrm{\left)}{\mathit{T}}_{\mathrm{rim}}\mathrm{\right(}\mathit{R}\mathrm{+}\mathit{dR}{\mathrm{)}}^{\mathrm{4}}\mathit{,}$(B.1)where ϵ(R) = (T_rim(R + dR)/T_rim(R))^β is the emissivity of the slab at R + dR, see in Appendix C. Note that here we neglect the irradiation of the neighboring slab at R + 2dR with a lower temperature than the slab at R + dR. To calculate T_rim(R), we first approximate T_rim(R + dR) with . This results in an ODE for T(R), which has the solution ${\mathit{T}}_{\mathrm{rim}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{rim}}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\exp \left[\mathrm{-}\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\beta }\mathrm{)}}\right)\mathrm{(}\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{0}}\mathrm{)}\right]\mathit{,}$(B.2)where R₀ is the radius of the disk inner edge. To take into consideration the additional irradiation of the rim at the opposite side, we assume that the temperature at the innermost slab is T_rim(R₀) = qT_int,irr(R₀), where q ≃ 1.2 according to results of our simulations done by 2D RADMC (Dullemond & Dominik 2004). Thus incorporating the additional heating by the disk rim, the disk interior temperature, previously given by Eq. (A.12), can be given by ${\mathit{T}}_{\mathrm{int}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\left(}{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{irr}}\mathrm{\right(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}\mathit{,}\mathrm{acc}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{rim}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{4}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}$(B.3)where the superposition rule is applied.

Appendix C: Dust emissivity

In this section the derivation of dust emissivity at a specific temperature is given, assuming optically thick environment. Let us define the dust emissivity at a temperature T_dust heated by a blackbody of the temperature T_irr as the ratio of the Rosseland mean opacities ${\mathit{ϵ}}_{\mathrm{dust}}\mathrm{=}\frac{{⟨\mathit{\kappa }\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}⟩}_{\mathrm{R}}}{{⟨\mathit{\kappa }\mathrm{\left(}{\mathit{T}}_{\mathrm{irr}}\mathrm{\right)}⟩}_{\mathrm{R}}}\mathrm{·}$(C.1)The Rosseland mean opacity is defined by ${⟨\mathit{\kappa }\mathrm{\left(}\mathit{T}\mathrm{\right)}⟩}_{\mathrm{R}}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{B}\mathrm{\left(}\mathit{\nu ,T}\mathrm{\right)}\mathit{/}\mathrm{d}\mathit{T}\mathrm{d}\mathit{\nu }}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{B}\mathrm{\left(}\mathit{\nu ,T}\mathrm{\right)}\mathit{/}\mathrm{d}\mathit{T\kappa }\mathrm{(}\mathit{\nu ,T}{\mathrm{)}}^{-1}\mathrm{d}\mathit{\nu }}\mathit{,}$(C.2)where B(ν,T) is the Planck function. The dust used in our model generates the κ(ν) = κ₀(ν/ν₀)^β opacity law, where β > 0. Note that wee use β = 1 throughout all our models (Rodmann et al. 2006). Following the calculations of Stahler & Palla (2005) in Appendix G, we assume that the dust absorption coefficient can be given by $\mathit{\kappa }\mathrm{\left(}\mathit{\nu ,T}\mathrm{\right)}\mathrm{=}{\mathit{\kappa }}_{\mathrm{0}}{\left(\frac{\mathit{kT}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathit{\beta }}{\mathit{x}}^{\mathit{\beta }}\mathit{,}$(C.3)where k and h are Boltzmann and Planck constants, respectively, and x = hν/kT. Supposing that $\mathit{\partial B}\mathrm{\left(}\mathit{\nu ,T}\mathrm{\right)}\mathit{/}\mathit{\partial T}\mathrm{d}\mathit{\nu }\mathrm{=}\mathit{Af}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x,}$(C.4)where A is an appropriate dimensional constant, and substituting this into the definition of Rosseland mean opacity, Eq. (C.2), leads to ${⟨\mathit{\kappa }\mathrm{\left(}\mathit{T}\mathrm{\right)}⟩}_{\mathrm{R}}\mathrm{=}{\mathit{\kappa }}_{\mathrm{0}}{\left(\frac{\mathit{kT}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathit{\beta }}\frac{\mathrm{\int }\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{\mathrm{\int }{\mathit{x}}^{\mathrm{-}\mathit{\beta }}\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}\mathrm{·}$(C.5)Because the quotient of integrals is a pure number, the Rosseland mean opacity is proportional to T^β. Applying (C.1) and (C.5) we find that the dust emissivity can be given by ${\mathit{ϵ}}_{\mathrm{dust}}\mathrm{=}{\left(\frac{{\mathit{T}}_{\mathrm{dust}}}{{\mathit{T}}_{\mathrm{irr}}}\right)}^{\mathit{\beta }}\mathrm{·}$(C.6)

Appendix D: Optical depth of the disk atmosphere

Assuming that the disk is observed with an inclination angle i, the monochromatic optical depth of the disk atmosphere is the sum of optical depths of dust and gas along the line of sight, i.e.: $\mathit{\tau }\mathrm{\left(}\mathit{\nu ,R,\phi }\mathit{,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\cos \mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\left(}{\mathit{\kappa }}_{\mathrm{d}}\mathrm{\right(}\mathit{\nu }\mathrm{\left)}{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\right(}\mathit{R}\mathrm{)}\mathrm{+}{\mathit{\kappa }}_{\mathrm{g}}\mathrm{(}\mathit{\nu ,R,\phi ,i}\mathrm{\left)}{\mathrm{\Sigma }}_{\mathit{g}}\mathrm{\right(}\mathit{R}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{,}$(D.1)where κ_d(ν) and κ_g(ν,R,φ,i) are the dust and gas opacity at frequency ν, while Σ_d(R) and Σ_g(R) are the dust and gas surface densities in the disk atmosphere. Note that because the dust opacity (κ_d) is taken to be uniform throughout the disk, the optical depth of the disk atmosphere depends on R and φ via the gas opacity characterized by the temperature distribution and the line-of-sight component of the orbital velocity of gas parcels. By definition the optical depth at the bottom of the disk atmosphere is unity at optical wavelengths along the stellar irradiation. Assuming that the mere opacity source at visual wavelengths is the dust, the dust surface-density of disk atmosphere is ${\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{\kappa }}_{\mathrm{V}}}\mathit{,}$(D.2)where κ_V is the overall opacity of dust at visual wavelengths. Assuming that the gas- and dust-mass ratio to the total mass (X_g and X_d, respectively) is constant throughout the disk (i.e. there are no vertical or radial variations in the mass ratios), the surface density of gas in the disk atmosphere is ${\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{X}}_{\mathrm{g}}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{X}}_{\mathrm{g}}}{{\mathit{X}}_{\mathrm{d}}}{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{,}$(D.3)using the dust surface-density Eq. (D.2). Thus, the optical depth in the disk atmosphere along the line of sight can be given by $\mathit{\tau }\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\cos \mathrm{\left(}\mathit{i}\mathrm{\right)}}\left({\mathit{\kappa }}_{\mathrm{d}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\frac{\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{\kappa }}_{\mathrm{V}}}\mathrm{+}\frac{{\mathit{X}}_{\mathrm{g}}}{{\mathit{X}}_{\mathrm{d}}}\frac{\mathit{\delta }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{\kappa }}_{\mathrm{V}}}{\mathit{\kappa }}_{\mathrm{g}}\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}\right)\mathit{,}$(D.4)where we used Eqs. (D.1)–(D.3).

Appendix E: Monochromatic opacities

In LTE the monochromatic opacity of the emitting gas at frequency ν with a molecular mass of m_CO due to transitions between states u → l can be given by $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{g}}\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{m}}_{\mathrm{CO}}\mathrm{8}\mathit{\pi }}\frac{\mathrm{1}}{\mathit{Q}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right(}\mathit{R}\mathrm{\left)}\mathrm{\right)}}{\left(\frac{\mathit{c}}{{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathrm{2}}{\mathit{A}}_{\mathrm{ul}}{\mathit{g}}_{\mathrm{u}}\\ & & \mathrm{\times }\left(\exp \left[\mathrm{-}\frac{{\mathit{E}}_{\mathrm{l}}}{\mathit{k}{\mathit{T}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right]\mathrm{-}\exp \left[\mathrm{-}\frac{{\mathit{E}}_{\mathrm{u}}}{\mathit{k}{\mathit{T}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right]\right)\\ & & \mathrm{\times }\mathrm{\Phi }\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}\end{array}$where ν₀ is the fundamental frequency of the transition, Q(T_atm(R)) is the partition sum at the gas temperature T_atm(R), A_ul and g_u are the probability of transition (i.e. the Einstein A coefficient of the given transition) and the statistical weight of the upper state, respectively, while c is the the light speed. The partition sum can be given by $\mathit{Q}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right(}\mathit{R}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}{\mathit{g}}_{\mathit{i}}\exp \left[\mathrm{-}\frac{{\mathit{E}}_{\mathit{i}}}{\mathit{k}{\mathit{T}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right]\mathit{,}$(E.3)where g_i and E_i are the statistical weight and energy level of the ith excitation state. In Eq. (E.1) the Φ(ν,R,φ,i) is the local intrinsic line profile originating by the natural thermal and the local turbulent broadening acting together. If the pressure of gas is negligible, the intrinsic line profile can be represented by a normalized Gauss function $\mathrm{\Phi }\mathrm{\left(}\mathit{\nu ,R,\phi ,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\sqrt{\mathit{\pi }}}\exp \left[\mathrm{-}{\left(\frac{\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}\mathit{R,\phi ,i}\mathrm{\right)}}{\mathit{\sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right)}^{\mathrm{2}}\right]\mathit{,}$(E.4)where σ(R) is the line width, Δν(R,φ,i) is the line center shift due to Doppler shift caused by the apparent motion of the gas parcels along the line of sight. The line width is determined by the natural thermal broadening ${\mathit{\sigma }}_{\mathrm{therm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\nu }}_{\mathrm{0}}}{\mathit{c}}\sqrt{\frac{\mathrm{2}\mathit{k}{\mathit{T}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{m}}_{\mathrm{CO}}}}\mathit{,}$(E.5)and the local turbulent broadening ${\mathit{\sigma }}_{\mathrm{turb}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\nu }}_{\mathrm{0}}}{\mathit{c}}\mathit{\chi }\sqrt{\frac{\mathit{\gamma k}{\mathit{T}}_{\mathrm{atm}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{m}}_{\mathrm{H}}}}\mathit{,}$(E.6)assuming that the speed of turbulent motions is χ times the local sound speed. In Eqs. (E.5) and (E.6) m_CO and m_H are the molecular masses of H and CO and γ is the adiabatic index of the main constituent of the gas, i.e. that of hydrogen. Considering the thermal and local turbulent broadening, the resulting profile will be the convolution of the two Gaussian line profiles, i.e. $\mathit{\sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathrm{therm}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{turb}}\mathrm{(}\mathit{R}{\mathrm{)}}^{\mathrm{2}}}\mathit{.}$(E.7)In a planet-free disk, in which the gas parcels are moving on circularly Keplerian orbits, the line center at ν₀ fundamental frequency shifts due to the Doppler shift is $\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}\mathit{R,\phi ,i}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\nu }}_{\mathrm{0}}}{\mathit{c}}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\mathit{R}}}\cos \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{.}$(E.8)Here we neglect that the massive planet and the host star are orbiting the common center of mass, instead the center of mass is set to the center of host star. Moreover, the influence of the gas pressure on the angular velocity of gas parcels, which causes slightly sub-Keplerian orbital velocities due to radial pressure support, is also not taken into account.

All Tables

All Figures