© ESO, 2014

1. Introduction

Observations of young stars of different ages reveal that protoplanetary disks evolve from optically thick, gas-rich disks to optically thin, gas-poor debris disks (e.g., see review by Williams & Cieza 2011). The transition between these two classes of objects is believed to occur relatively fast (10⁵ years) compared to the disk’s life time (5−10 Myr) (e.g., Cieza et al. 2007; Damjanov et al. 2007; Currie & Sicilia-Aguilar 2011). Infrared space observatories have unveiled several young stars with spectral energy distributions (SED) characterized by IR-excess at >10 μm and significantly reduced excesses at shorter wavelengths (e.g., Strom et al. 1989; Calvet et al. 2005; Sicilia-Aguilar et al. 2006; Hernández et al. 2006; Espaillat et al. 2007; Fang et al. 2009; Merín et al. 2010; Rebull et al. 2010; Cieza et al. 2010). These sources are commonly named “transition disks”, because they are believed to be the objects in the transition phase between a young star with an optically thick gas-rich disk and a star with an optically thin gas-poor debris disk.

Generally, the lack of strong emission in the near or mid-IR SED is interpreted as evidence of a gap or cavity from a few up to tens of AU in the disk. Follow-up imaging of several transition disks with sub-mm interferometers have confirmed that transition disks indeed have a deficit of sub-mm continuum emission at a few mJy levels inside tens of AU (e.g., Piétu et al. 2007; Brown et al. 2009; Hughes et al. 2009; Andrews et al. 2011; Isella et al. 2012), thus providing further evidence for such inner disk cavities, at least on the large dust grains.

Several scenarios have been discussed in the literature to explain the presence of a cavity and the transitional disk SED shape: grain growth (e.g., Dullemond & Dominik 2005; Birnstiel et al. 2012), migrating giant planets that opened a gap (e.g., Varnière et al. 2006; Zhu et al. 2011; Dodson-Robinson & Salyk 2011), dust filtration by a giant planet (e.g., Rice et al. 2006; Zhu et al. 2012; Pinilla et al. 2012), disk dissipation due to a photoevaporative disk wind (e.g., Alexander & Armitage 2007; Owen et al. 2012), and dust free inner holes due to radiation pressure (Chiang & Murray-Clay 2007). One important step toward understanding the transitional disk phenomenon is to understand how the gas and dust structures compare. Do they follow each other in density structure? Do they thermally de-couple? Does the gas-to-dust ratio stay constant as a function of the radius?

Several young stars with transition disks are emission-line stars that exhibit signs of accretion (e.g., Najita et al. 2007b; Muzerolle et al. 2010; Cieza et al. 2012). As a consequence, the cavities imaged in the sub-mm should have gas as the material from the optically thick outer disk flows through the cavity to be accreted by the central object (e.g., Lubow & D’Angelo 2006). In fact, several transition disks display emission of CO at 4.7μm (e.g., Salyk et al. 2009) a common tracer of warm gas in the inner disk (e.g., Najita et al. 2007a). Furthermore, in a few transition disks, CO ro-vibrational emission has been spatially resolved up to distances of tens of AU (Pontoppidan et al. 2008). Moreover, recent studies of dust-scattered light in the near-IR have detected emission from small dust grains inside the cavities that were previously imaged with sub-mm interferometry, with the additional and remarkable property that no sharp edge is seen at the location of the sub-mm cavity’s inner radius (e.g. Muto et al. 2012; Garufi et al. 2013). This suggests that what is observed is likely a change in the dust size distribution inside the cavity (Dong et al. 2012). Near-IR interferometry provides evidence of dust emission (and inhomogeneities) inside the dust cavity of transition disks (e.g., Kraus et al. 2013). Furthermore, recent ALMA observations have spatially resolved gas inside the cavity of transition disks (Casassus et al. 2013; Bruderer et al. 2013).

In summary, the observational evidence indicating that the dust cavities of transitions disks are indeed not empty has been steadily growing during past years. This opens interesting questions such as how much gas and dust are present in the cavity; what is their distribution and chemical composition; if dust is present, what is its size distribution; if there is a gap in the gas or in the dust, how large is it?

The goal of this paper is to address the problem of the gas and dust structures in transition disks by performing a detailed study of the transition disk HD 135344B (SAO 206462). Our aim is to constrain the disk structure by modeling simultaneously and in a coherent way multiwavelength & multi-instrument gas and dust observations of the inner and the outer disk. Most specifically, we investigate gas and dust disk content and structure inside the cavity by simultaneously modeling the SED, new VLTI-PIONIER near-IR interferometry data, the CO 4.7 μm emission, and constraints from sub-mm continuum interferometry. Furthermore, we seek to constrain the gas mass and the gas-to-dust ratio in the outer disk by employing the SED, the [O i] 63 μm and 145 μm lines, and (sub-)mm observations.

HD 135344B is an accreting (2 × 10^-8 M_⊙/yr, Sitko et al. 2012) F4V young star (8$\begin{array}{c}\mathrm{+}\mathrm{8}\\ -4\end{array}$ Myr, van Boekel et al. 2005) that has a “pre-transition disk” (i.e., a transition disk with near-IR excess). It has the remarkable characteristic of exhibiting CO 4.7 μm emission extending tens of AU inside the 45 AU sub-mm dust cavity (Pontoppidan et al. 2008) and scattered light emission from small dust grains down to 28 AU (Muto et al. 2012; Garufi et al. 2013). HD 135344B display spiral structures (Muto et al. 2012; Garufi et al. 2013) and asymmetries in its outer disk dust emission (Brown et al. 2009; Andrews et al. 2011; Pérez et al. 2014). It is a transition disk that exhibits variability in its near-IR SED and in optical and near-IR line emission on time scales of months (Sitko et al. 2012). It is a source on which [O i] 63 μm emission has been detected in Herschel/PACS observations (Meeus et al. 2012). Finally, it is a transition disk without close-in low-mass stellar companions (Vicente et al. 2011).

We start the paper by a brief summary of current observational constraints on the disk of HD 135344B in Sect. 2. Then, in Sect. 3, we describe the modeling tools and the general modeling procedure. In Sect. 4, we present and discuss the different models that were tested. In Sect. 5, we examine the disk constraints derived from the final disk model. In Sect. 6, we discuss our results in the context of the study of transitional disks and planet formation. Finally a summary and conclusions are presented in Sect. 7.

2. Observational constraints

HD 135344B has been observed with a diversity of instruments and techniques from UV to mm wavelengths. In Table 1, we present a concise summary of relevant previous observational constraints on HD 135344B. The most important observational constraints for this study are:

• the SED has a dip in the 10−50μm region (Brown et al. 2007);

• the 10 μm silicate feature is absent (Geers et al. 2006);

• the inner radius of the sub-mm cavity is 46±5 AU (Brown et al. 2009; Andrews et al. 2011);

• the 3σ continuum flux upper limit at 870 μm is 10.5 mJy for a 0.24′′ × 0.5′′ beam centered on the star (Andrews et al. 2011);

• CO ro-vibrational emission extends at least to 15 AU inside the sub-mm dust cavity (25 AU for d = 140 pc) (Pontoppidan et al. 2008)¹;

• near-IR dust scattering images reveal small dust inside the sub-mm cavity down to 28 AU with a smooth surface brightness and no discontinuity at 45 AU (Muto et al. 2012; Garufi et al. 2013).

3. Modeling

3.1. Stellar parameters

We provide a detailed discussion of the derivation of the stellar parameters in the Appendix A.2. For our models we used a star with T_eff = 6620 K (F4V), a stellar radius of 2.1 R_⊙, a mass of 1.65 M_⊙, and A_v = 0.4. Owing to uncertainties in the photometry, the spectral type and, especially, the distance (±20 pc) the stellar radius can vary ±0.3 R_⊙, and the stellar mass does so ±0.1 M_⊙. The stellar UV spectrum² was parametrized with a fractional UV excess f_UV = L_UV/L_⋆ equal to 0.001 (L_UV is defined as the luminosity between 91.2 and 250 nm) and F_ν ∝ ν^γ with γ = −2.15. This parametrized stellar UV spectrum is included in the ProDiMo gas heating calculation in addition to the interstellar UV field.

3.2. Disk inclination and position angle

Several data sets give different inclinations for HD 135344B: CO 4.7 μm emission suggests i = 14° ± 4° (Pontoppidan et al. 2008, for M = 1.65 M_⊙); the CO J = 3−2 line at 870 μm indicates i = 11° ± 2° (Dent et al. 2005); sub-mm continuum imaging points toward i = 12° (Andrews et al. 2011) and 21° (Brown et al. 2009); mm interferometry suggests i = 11° (Lyo et al. 2011); near-IR imaging sets an upper limit of i = 20°Grady et al. (2009); mid-IR imaging suggests i = 46° ± 5° (Doucet et al. 2006) and 45°± 5° (Mariñas et al. 2011). For our analysis of the SED and line profiles, we used i = 14°. We used this value because it is derived from a simultaneous fit to the spectro-astrometry signature of the CO 4.7 μm line in three-slit position angles and the spectrally resolved CO 4.7 μm line profile, and because it is consistent with the inclination derived from the observations of CO gas in the outer disk³.

Several estimations exist for the disk’s position angle (PA). Pontoppidan et al. (2008) suggest a PA = 56°± 2° based on spectro-astrometry of the CO 4.7 μm emission; Grady et al. (2009) suggest a PA = 55°± 5° based on near-IR scattered light imaging; Andrews et al. (2011), Brown et al. (2009), and Lyo et al. (2011) based on (sub-)mm interferometry continuum observations suggest a PA = 64°, 55°, and 64° respectively. As we aim to compare our model with the CO 4.7 μm lines we employed a PA of 56°.

3.3. Description of the general modeling procedure

We aimed to

• 1.

  simultaneously fit the SED, the line profile of the CO ν = 1−0 P(10) line at 4.7545 μm, and the near-IR PIONIER visibilities and closure phases;

• 2.

  reproduce within a factor of a few the detected line fluxes of CO P(10), [O i] at 63 μm, ¹²CO J = 3−2 at 870 μm, and ¹²CO J = 2−1 at 1.27 mm;

• 3.

  obtain line fluxes below the upper limits for [O i] at 145 μm, [C ii] at 157 μm, H₂1−0 S(1) at 2.12 μm, H₂0−0 S(1) at 17 μm, and CO and H₂O infrared lines covered by Herschel (see Tables A.2 and A.3);

• 4.

  obtain a 870 μm flux inside a beam of 4′′×0.5′′(33×70 AU) centered on the star lower than 10.5 mJy.

To be consistent with the sub-mm continuum constraints, the inner radius of the outer disk’s large and small grains was set to 45 AU. Later in the modeling process, the inner radius of the outer disk’s small grains was allowed to extend down to 30 AU to be consistent with the near-IR scattered-light constraints.

In Fig. 1, we provide an schematic overview of the modeling procedure. The modeling starts by assuming a dust composition, a dust size distribution, a dust mass, and a gas-to-dust ratio, for the inner and the outer disk using a parametric disk (see details in next section). Then the Monte Carlo dust radiative transfer code MCFOST (Pinte et al. 2006, 2009) is used to calculate the dust’s thermal and density structure, T_dust(r,z), n_dust(r,z), and the mean radiation field J_ν(r,z). MCFOST is employed to compute a synthetic SED, and continuum images in the H band and at 870 μm for comparison with VLTI/PIONIER and SMA continuum observations.

The MCFOST disk structure and radiation field is then used as input for the thermochemical radiative transfer code ProDiMo (Woitke et al. 2009). ProDiMo calculates the disk chemistry, the gas and dust heating and cooling, and the energy level populations of molecules and atoms. ProDiMo computes the synthetic line fluxes and channel maps for gas emission lines from the optical to the mm using escape probability.

Details of the coupling between the codes MCFOST and ProDiMo and brief summaries of technical details of both codes are provided in Woitke et al. (2010), Pinte et al. (2010), and Kamp et al. (2011). Details of the implementation of CO ro-vibrational emission within ProDiMo are given in Thi et al. (2013).

The modeling was performed without computing the hydrostatic equilibrium. The disk gas and dust density distribution (thus the gas-to-dust ratio) was set as input. We chose this approach to explore a large fraction of the parameter space. Computing the hydrostatic equilibrium would have resulted in a running time an order of magnitude longer and would have prevented us from performing a large exploration of the parameter space. We have verified that our last model, Model 5, is consistent with hydrostatic equilibrium; i.e., the scale height we used as input is consistent with ${\mathit{h}}_{\mathrm{hydro}}\mathrm{=}\sqrt{\mathit{k}{\mathit{T}}_{\mathrm{gas}\hspace{0.17em}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}}\hspace{0.17em}{\mathit{r}}^{\mathrm{3}}\mathit{/}\mathit{GM\mu }}$.

We started our exploration with simple models, and complexity in the model was only added when we were not able to fit some of the observations. The model was refined sequentially following this protocol:

• 1.

  fit the SED,

• 2.

  fit the CO P(10) ro-vibrational line-profile,

• 3.

  check consistency with SMA 870 μm continuum upper limits inside the cavity,

• 4.

  describe near-IR visibilities,

• 5.

  produce an [O i] 63 μm line flux consistent with the Herschel observations,

• 6.

  refine the model to account for the latest near-IR polarization images (Garufi et al. 2013).

For each family of models, a large portion of the parameter space was explored by a combination of hand exploration, grid modeling, and simplex optimization algorithms. In total we tested ~60 000 MCFOST models and ~4000 ProDiMo models.

3.4. Disk structure and MCFOST parameters

The philosophy of MCFOST is to model the disk as composed of individual zones that may overlap. Each zone is a disk extending from an inner radius R_in to an outer radius R_out, with a scale height h₀ at a reference radius r₀, flaring exponent (β), and surface density exponent (q). Each zone has a Gaussian vertical density profile ρ(r,z) = ρ₀(r) exp(−z²/2h(r)²), a power-law surface density Σ(r) = Σ₀(r/r₀)^q, and a scale height h(r) = h₀(r/r₀)^β, where β is the flaring exponent, r the radial coordinate in the equatorial plane, and h₀ the scale height at the reference radius r₀.

For each zone, the dust composition, the dust size distribution, the dust mass, and the gas-to-dust ratio are defined independently. One zone can include several dust components, such as silicate grains, carbonaceous grains, or PAHs. Zones can have Gaussian inner edges of size five times the parameter “edge”. In our models, for each zone, the gas-to-dust ratio is constant in the vertical direction.

Dust grains are defined as homogeneous and spherical particles with sizes distributed according to the power law dn(a) ∝ a^p da, with a_min and a_max the minimum and maximum sizes of the grains. We used the standard value p = −3.5. Extinction and scattering opacities, scattering phase functions, and Mueller matrices are calculated using the Mie theory. We employed the astronomical silicates opacities of Draine & Lee (1984), the neutral PAH opacities computed by B. T. Draine (compiled from Laor & Draine 1993; Draine & Lee 1984; Li & Draine 2001), and the amorphous carbonaceous dust optical constants derived by Li & Greenberg (1997).

In summary, the free parameters of a model are the number of zones used to describe a disk, and for each zone R_in, R_out, p, β, q, h/r, edge, M_dust, gas-to-dust ratio, and dust composition, and for each dust species a_min, and a_max.

Sub-mm continuum images were convolved with a 0.24″ × 0.5″ beam for comparison with SMA observations.

3.5. ProDiMo parameters

We used the cosmic ray ionization rate of H₂ of 5 × 10^-17 s^-1. The incident vertical UV was set to the interstellar medium value. Non-thermal broadening was set to 0.15 km s^-1. We used the UMIST chemical network (9 elements, 71 species connected through 950 reactions neutral-neutral, ion-molecule, photoreactions, cosmic ray reactions, and absorption & desorption of CO, CO₂, H₂O, NH₃, CH₄, see Woitke et al. 2009). UV fluorescence is included for [O i] and [C ii] but not for CO. UV CO fluorescence has been included as a test later in the paper. A list of the species we used is provided in Woitke et al. (2009).

3.6. CO 4.7 μm data and slit effects

For our analysis, we downloaded the 1D CO 4.7 μm reduced spectrum used in Pontoppidan et al. (2008)⁴. We flux-calibrated the spectrum using the Spitzer flux at 4.7 μm (see Appendix A for the photometry) and measured the CO line fluxes by fitting a Gaussian to their line profile. The CO P(1) to P(11) line profiles are available within the CRIRES spectrum. After subtracting the continuum and normalization by the peak flux, the lines have the same line profile within the errors (see Fig. A.3). We selected the CO P(10) line for detailed modeling because the line profile is complete, is weakly affected by nearby strong telluric absorption lines, and has good signal-to-noise ratio (S/N). The line CO P(10) line flux is 1.5 × 10^-17 W m^-2 with an error on the order of 20% owing to slit losses and systematics.

To compare the ProDiMo predictions with the observed CRIRES CO 4.7 μm spectra, the effects of the slit width and orientation, the observing conditions, and spectral resolution needed to be taken into account.

ProDiMo generates channel maps data cubes; i.e., for each velocity bin, a 2D image of CO emission is generated (see details in Hein Bertelsen et al. 2014). The pixel size of the CRIRES detector in the spatial direction is 0.086′′, at a distance of 140 pc that corresponds to a pixel size of ~12 AU. We generated the ProDiMo data cubes with a pixel size of 2 AU.

The ProDiMo data cubes were first convolved in the spatial direction with a Gaussian PSF of FWHM = 180 mas. Then they were convolved with a Gaussian of FWHM = 3.3 km s^-1 in the wavelength direction to simulate R = 90 000 of the observations. A mask of 0.2′′ (28 AU) with the correct PA was applied in each channel to mimic the slit. The fluxes were added in the direction perpendicular to the slit to generate a 2D spectrum.

In the 2D spectrum the centroid of the photocenter at line velocity (υ) was calculated using (Pontoppidan et al. 2008) $\mathit{X}\mathrm{\left(}\mathit{\upsilon }\mathrm{\right)}\mathrm{=}\mathit{K}\frac{{\sum }_{\mathit{i}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right(}\mathit{\upsilon }\mathrm{)}\mathrm{-}{\mathit{x}}_{\mathrm{0}}\mathrm{)}{\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{\upsilon }\mathrm{\right)}}{{\sum }_{\mathit{i}}{\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{\upsilon }\mathrm{\right)}}\mathrm{\left(}\mathrm{pixels}\mathrm{\right)}$(1)where x_i − x₀ is the center of pixel i relative the continuum centroid position, F_i(υ) is the flux on the pixel i, and K is a correction factor to take into account that not all of the source flux is enclosed in the aperture (Pontoppidan et al. 2008). We employed K = 1.3 used by Pontoppidan et al. (2008) in their analysis of the spectro-astrometry signal of HD 135344B to be able to compare our models to their measurements. A 1D spectrum is further obtained by summing the pixels in the spatial direction. To compare the observed and synthetic line profiles, the spectra were normalized by dividing them by the median of the continuum, then after continuum subtraction, the profile was renormalized by dividing it by its maximum flux. In this way the continuum is always at 0 and the line peak is always at 1. The integrated line fluxes of the model (taking the slit effects into account) and the line profiles were compared to the observations separately.

4. Modeling results

Before presenting the family of models 5, the model that best describes the observations, we briefly describe the sequence of families of models we tested to illustrate the reasoning that lead us to the final model. A family of models is defined primarily by its dust properties (i.e., composition and location). We have extensively explored the parameter space in each family of models. However, we limit the discussion to the parameter space on the family of models 5. In Table 4 we present the details of a representative example of each family of models. The different families of models show us the close relation that exists between the properties assumed for the dust (size distribution, composition, and mass) at the beginning of the modeling procedure and the gas lines obtained.

4.1. Family of models 1: a disk only composed of astronomical silicate grains

As starting point we assumed the simplest model, which has an inner and outer disk composed of 100% astronomical silicates and a gas-to-dust ratio of 100 in the whole disk. The solutions all converged to the same disk structure, namely a narrow ring of dust at 0.16 to 0.21 AU with a dust mass of a few 10^-10 M_⊙ with grains larger than 10 μm (to reproduce the lack of 10 μm silicate feature), followed by a large gap of 45 AU and an outer disk from 45 to 200 AU with dust mass of 10^-4 M_⊙. By changing the disk’s geometrical parameters, we found a disk able to reproduce most of the observed gas line fluxes (see Model 1a in Tables 4 and 5). However, because in this model the warmgas is located in the narrow innermost region, the CO ro-vibrational emission produced has a broad double-peaked profile with FWHM ~ 50 km s^-1 that is inconsistent with the observed FWHM ~ 15 km s^-1 (see Model 1a in Fig. 3).

With the aim of producing CO ro-vibrational emission extending tens of AU and a narrower line profile, in Model 1b, we introduced gas (and dust) between 0.21 and 45 AU. We found that although up to 10^-8 M_⊙ of dust and 10^-6 M_⊙ of gas can be introduced inside the cavity and still fit the SED, the CO ro-vibrational line always displayed a broad double peaked profile with a FHWM of ~50 km s^-1. We tested many flat and flared disks filling the cavity (see Model 1b in Fig. 3), but the same result was found in all cases. The optically thick innermost disk of silicates required to fit the near-IR SED shields the material inside the cavity and does not permit the gas at several AU to be warm enough to contribute significantly to the CO ro-vibrational lines. A change in the dust’s optical properties of the inner disk is required to enable ro-vibrational CO emission at tens of AU.

4.1.1. PAH content

PAHs are an important ingredient for calculating the gas heating. To constrain the PAHs, we employed the Spitzer/IRS spectrum⁵. PAHs are implemented as a second dust component in MCFOST. Their properties (abundance and size) are passed to ProDiMo in order to compute the gas heating due to the photoelectric effect.

Initially, we found that the observed 11.2 μm PAH feature could be reproduced with neutral PAHs with 21 carbon atoms (NC = 21) and 0.01% of the dust in the outer disk in the form of PAHs. However, such PAHs generated too much gas heating that translated into [O i] 63 μm line fluxes that are three to ten times stronger than the observations for almost all the models reproducing the SED. To have an [O i] 63 μm emission compatible within a factor 3 of the observations, we increased the PAH size to NC = 100 and used a 0.01% fraction in mass of the dust in the form of PAHs (f_PAH = 6.0 × 10^-3 for a gas-to-dust ratio of 100, see Fig. 2).

4.2. Family of models 2: an inner disk composed of a uniform mixture of amorphous carbon and astronomical silicates

To modify the continuum optical depth of the inner disk and to make it possible for the gas at larger radii (R ≫ 1 AU) to contribute to the CO 4.7 μm emission, we introduced amorphous carbonaceous grains in the dust mixture in the family of models 2. Amorphous carbonaceous grains are commonly employed to fit SEDs of T Tauri and Herbig Ae/Be stars. We assumed first a carbonaceous/silicate grains ratio constant with radius. The fraction of carbon grains in the dust mixture sets the extension of the inner disk that fits the near-IR SED. The greater the fraction of carbon grains, the larger the inner disk that reproduces the near-IR continuum. With a 25% fraction of carbonaceous grains, an inner disk extending from 0.18 to 20 AU fits the near-IR SED (see Model 2a in Table 4). With a 5% carbon fraction, an inner disk extending from 0.18 to 3 AU is required (Model 2b).

The introduction of amorphous carbon grains reduced the total dust mass required to fit the near-IR SED by a factor of a few, allowed for a dust size distribution with smaller grains in the inner disk (and fit the lack of 10 μm silicate feature), and drastically changed the inner disk’s optical continuum depth. These changes resulted in an optically thin inner disk at 4.7 μm, thus a CO ro-vibrational emission entirely dominated by emission from the inner rim of the outer disk, and a very narrow single peaked CO ro-vibrational line profile with a width of a few km s^-1, inconsistent with the observations (see Fig. 3).

4.3. Family of models 3: an inner disk with a radial-dependent mixture of carbon and silicates

The results of the families of models 1 and 2 suggested that the solution was an intermediate inner disk structure between a narrow disk of 100% silicates and an extended inner disk with a large amount of carbonaceous grains. In the family of models 3, we allowed for a radial-dependent carbonaceous/silicate grains ratio.

We found that an inner disk of tens of AU composed of silicates, but enriched with a small amount (10^-12 M_⊙) of amorphous carbon grains at small radii (0.2 <R< 0.26 AU) is able to reproduce the near-IR SED, while at the same time it makes it possible for the CO located at several AU to contribute significantly to the CO ro-vibrational line flux, thus reproducing the observed CO P(10) line profile (see Fig. 3).

Assuming a gas-to-dust ratio of 100 for the whole disk and a carbon-enriched inner disk, we tested a large number of models (~50 000 MCFOST, ~1000 ProDiMo) by varying the geometry (i.e., H/R, flaring, surface density exponent) and dust mass for the inner and outer disk. We found that the CO ro-vibrational line profile is reproduced by an inner disk extending tens of AU in which the surface density is flat or increases as a function of the radius (i.e., a power-law surface density with a positive exponent). Solutions with flared and anti-flared outer disks were found.

Several Models 3 reproduced the CO P(10) flux and line profile, the H₂ infrared lines upper-limits, and the CO sub-mm and mm line fluxes (see one example in Tables 4 and 5). However, in all Models 3 calculated (covering a wide range of geometries), the emission of [OI] 63 and 145 μm, and the [C ii] line at 157 μm were overpredicted by factors over five. We tested several options available within ProDiMo (see Fig. 4), such as global gas-to-dust ratios lower than 100, H₂O ro-vibrational cooling (pure H₂O rotational cooling is calculated by default), lower metallicity, no PAHs, extremely low O/C abundances, but in most of the cases the [O i] 63 and 145 μm lines were still too strong. We found that the only effective ways to significantly decrease the [O i] 63 and 145 μm line fluxes were to assume T_gas = T_dust or to significantly decrease the gas-to-dust ratio. The first is an extreme case of gas cooling that is unrealistic, as we know that T_gas>T_dust in the disk’s surface layer (e.g., Kamp & Dullemond 2004). The second had the limitation that when the global gas-to-dust ratio was low enough to describe the [O i] 63 μm line, the CO ro-vibrational line was two orders of magnitude weaker. This last result indicated that the gas-to-dust ratio should be different for the inner and outer disk.

Another limitation of the family of models 3 was the fit to the PIONIER near-IR visibilities. In a forthcoming paper (Benisty et al., in prep.), ESO-VLTI PIONIER near-IR observations of HD 135344B and other transition disks will be presented in detail. In Fig. 5, we present the observed VLTI/PIONIER squared visibilities at 1.6 μm and the predicted visibilities of Model 3. We see clearly that an inner disk starting at R_in = 0.2 AU is not compatible with the interferometry data. The near-IR continuum in Model 3 is produced too far out.

4.4. Family of models 4: Carbon grains inside the silicate sublimation radius, different gas-to-dust ratios for the inner and outer disk, and dust settling in the outer disk

4.4.1. The location of the carbonaceous dust

In the family of models 3, the carbonaceous grains were located at R> 0.2 AU to be consistent with the assumption that the temperature in the inner rim should not be higher than the silicate grains sublimation temperature (T ~ 1500 K, for n_H ≳ 10¹⁶ cm^-3, e.g., Helling et al. 2001). However, carbonaceous grains can survive temperatures higher than 1500 K and up to 2000 K (see for example, Kobayashi et al. 2011). Therefore, one interesting possibility for physically justifying the inner-disk carbon enrichment is that the carbonaceous grains are located inside 0.2 AU, the region where silicate grains sublimate. The survival of carbonaceous grains in this region is discussed in Sect. 6.2.

Therefore for the family of models 4, we located the carbonaceous grain component at 0.08 <R< 0.2 AU, R_in equal the corotation radius, and R_out equals the silicates sublimation radius. The astronomical silicate-grain component of the inner disk was set to start at 0.2 AU. With this inner disk configuration we obtained an excellent match to SED and the PIONIER visibilities and closure phases (see Fig. 5). In our models the temperature at the innermost radius is ~2000 K.

4.4.2. Decreasing the [O I] 63 μm emission

One limitation of the family of models 3 was that it produced [O i] 63 μm line fluxes that are too strong. The [O i] 63 μm line is emitted partly by gas inside the cavity between 10 and 45 AU, but principally by gas in the outer disk between 45 and 60 AU. Since the surface density and gas mass in the inner disk are set by the CO ro-vibrational line, to lower the [O i] line fluxes, three modifications to the family of models 3 were introduced in the family of models 4:

• 1.

  the gas-to-dust ratio in the outer disk was allowed to be lower than 100,

• 2.

  the surface density of the outer disk needed to be shallower, (We obtained good solutions with q =− 1.0, similar to the surface density power law exponent found by Andrews et al. (2011) from SMA sub-mm continuum observations).

• 3.

  the scale height of the outer disk needed to be lower than 10%.

These changes in the outer disk geometry led to a bad fit of SED at λ> 15 μm. To resolve this, a fourth modification was introduced in the outer disk by splitting it into two superposed disks. A first disk with lower H/R (5%), with 80% of the gas and dust mass, and with a dust population of large grains (0.01 <a< 1000 μm); and a second disk with higher H/R (0.8−0.13), with the 20% remaining gas and dust mass, and with a dust population of smaller grains (0.01 <a< 10 μm). This modification aims to keep most of the gas mass of the outer disk at low H/R (to fit the [O i] 63 μm line), while allowing some small dust particles to be present in an extended outer disk atmosphere at higher H/R to fit the SED at λ> 15 μm. This two-layered outer disk echoes the expected effect by dust coagulation and sedimentation (i.e., large grains closer to the midplane).

Finally, for the inner disk, to address the fact that scattered light imaging (Muto et al. 2012; Garufi et al. 2013) revealed material inside the cavity down to 28 AU, we set, as a first approximation, the inner disk to have an outer radius of 45 AU, such that the whole cavity is replenished with gas, and we allowed the gas-to-dust ratio to be over 100 in the inner disk to be able to have sufficiently high CO 4.7 μm emission.

With this ensemble of modifications, we found a family of disk models able to simultaneously describe the SED, the CO P(10) line profile, the line fluxes of [O i] 63 μm (within a factor 2), CO P(10) (within a factor 4), ¹²CO 3−2, ¹²CO 2−1, and the upper limits of [O i] at 145 μm, [C ii] at 157 μm, and the H₂ lines in the near and mid-IR. The properties of a representative model of the family of models 4 is described in Table 4, and the line fluxes predicted are presented in Table 5.

The family of models 4 confirmed the result that to describe the CO P(10) line profile, the gas in the inner disk should be distributed with a surface density increasing as a function of the radius, and indicated that to describe the [O i] 63 μm line flux, the gas-to-dust ratio in the outer disk should be much lower than 100. The best match to the [O i] 63 μm line flux was obtained by a gas-to-ratio below 10; nevertheless, gas-to-dust ratios up to 40 provided [O i] 63 μm fluxes within a factor 2 of the observations.

4.5. Family of models 5: introducing recent constraints from polarized and mid-IR imaging

In a recent work Garufi et al. (2013) has shown that polarized scattered light in HD 135344B drops significantly at a radius of 28 AU. A lack of polarized emission can be the signature of either a lack of material or a change in the illumination of material inside the cavity. Similarly, Maaskant et al. (2013) propose a cavity of size 30 AU based on modeling the SED and mid-IR imaging. The detection of scattered light down to 28 AU, inside the sub-mm dust cavity of 45 AU, indicates a different spatial location for the small and large dust grains in the outer disk.

To account for the recent results from Garufi et al. (2013), we slightly modified the disk structure of the family of models 4. First, we shortened the outer radius of the inner disk to 30 AU and decreased its dust mass to keep it consistent with the SMA 870 μm upper limit. Second, we extended the small particle component of the outer disk down to 30 AU. The large grain component of the outer disk was kept at 45 AU to account for the sub-mm 870 μm imaging constraints.

With these modifications, we obtained a disk structure compatible with the Garufi disk structure, while keeping the fit to the SED and the gas lines. Our model has gas and (some) dust inside 30 AU to account for the extended CO ro-vibrational line observed. However, the amount of dust at R< 30 AU (<10^-7 M_⊙) is much less than the amount of dust in the outer disk (2 × 10^-4 M_⊙). The properties of a representative model of the family of models 5 is described in Table 4, and the line fluxes predicted are presented in Table 5. In Fig. 6, we present the synthetic SED, the predicted CO P(10) line profile and spectro-astrometry signal, and the plots describing the optical depth of the line and of the continuum, the cumulated line flux intensity, and the number density of the species as a function of the radius for the CO P(10), [O i] 63 μm, and CO 3-2 lines at 870 μm lines. Figure 6 also presents a plot of a surface density of the gas and the dust as a function of the radius. In Fig. A.4, we present the plots of the number density of H as a function of the radius, the dust and gas temperature in the disk, and similar plots to those in Fig. 6 but for the [O i] 145 μm, [C ii] 157 μm, and H₂ 0-0 S(1) 17 μm lines.

5. Disk structure constraints derived from Model 5

The interest of performing multi-instrument modeling is to use the constraints obtained from different gas and dust tracers to break the model degeneracies and narrow down the parameter space of possible solutions. Ideally, a Bayesian analysis of a large number of models covering a significant fraction of the parameter space should be performed. However, owing to the prohibitive amount of computing time that this kind of analysis would require when the heating and cooling balance and the chemistry calculation are included, we limit our discussion to the parameter space surrounding Model 5, the best solution found. Our solution is a model that reproduces most of the constraints imposed by observations; however, there is the possibility that the solution is not unique.

5.1. Compromises during the modeling procedure

Ideally, a model should be able to reproduce all the observations available. However, to be able to reproduce most of the observations simultaneously, we needed to relax the perfect fit for a few of them. ProDiMo includes a larger number of physical processes, but, not all the physics are included in the code. Furthermore, our models are axisymmetric, HD 135344B is known to display a spiral structure (Muto et al. 2012; Garufi et al. 2013). Not including these spiral structures may have an impact, or not, on the integrated line fluxes that we are fitting. But, we are modeling integrated quantities and spiral arms produce local changes, so, it is not clear whether their impact would be dramatic.

The first compromise is the fit to near-IR SED at 8−10 μm. Our model is fainter than the observations at those wavelengths. For an innermost carbon disk that reproduces the near-IR interferometry data and the SED near-IR excess, one can construct a disk between 0.2 to 30 AU with sufficient dust mass for a given silicate dust size distribution, such that the 8 to 20 μm excess is well reproduced (see for example Model 3). However, the upper limit of the continuum emission at 870 μm inside 30 AU sets a stringent limit on the amount of dust mass that can be put in the inner disk. Higher masses can be achieved using smaller grain size distributions; however, when the dust size distribution is dominated by dust with a_min< 10 μm, then the silicate feature appears. A solution for better fitting the near-IR SED might be to introduce an additional zone at higher H/R in the inner disk with small grains (a_min< 10 μm) with 1% to 10% of the inner disk’s dust mass.

The second compromise is the absolute flux of the CO ro-vibrational lines. A cavity without discontinuity in the gas implies that around half of the [O i] 63 μm and 145 μm lines will be emitted inside the cavity (see Fig. 6). Therefore, when the gas mass or temperature is increased to better fit the CO P(10) flux, the flux of the [O i] lines will also increase. Since the [O i] 63 and 145 μm excitation is better understood and tested than CO ro-vibrational excitation (the collision rates are known only within a factor ten (Thi et al. 2013)), and because the [O i] 63 μm line can be used to constrain the gas mass, we gave priority to reproduce the [O i] 63 μm and 145 μm line fluxes and upper limits simultaneously with the CO P(10) line profile over fitting the CO P(10) integrated line flux. Our models underpredict the CO P(10) line flux, hence the spectro-astrometry signature.

5.2. Inner disk surface density exponent

The CO P(10) line profile strongly depends on the surface density power law exponent of the inner disk gas. In Fig. 7, we display the calculated CO P(10) line profile after taking the slit effects into account for disks with inner and outer disk H/R ranging from 0.07 to 0.15, and q ranging from −0.4 to +0.6. The CO P(10) line profile is best described by slightly positive power-law surface density exponents, in other words, by a surface density that increases as a function of the radius. Power law exponents q smaller than 0.0 produce line profiles that are too broad, and q larger than +0.4 produce line profiles that are too narrow. Keeping in mind the uncertainties and simplifications in the ProDiMo code, we can exclude gas surface density profiles with negative exponents for the inner disk. A steep surface density profile with a negative exponent will assign too much gas to the inner few AU of the disk. This would result in too much hot CO gas close to the star emitting at high-velocities; hence, a broad line-profile (see upper row in Fig. 7) that is inconsistent with the line profile and spectro-astrometric signal observed, which shows that the CO ro-vibrational emission comes from radii greater than a few AU.

This result is relatively independent of the scale height assumed for the inner disk. Here we have assumed that the CO 4.7 μm line is produced entirely by disk emission. Although a contribution from a weak disk wind cannot be ruled out (the object displays asymmetric oxygen emission at 6300 Å, van der Plas et al. 2008), the symmetric spectro-astrometry signature detected by Pontoppidan et al. (2008) at three slit positions favors a dominant contribution from disk emission.

Even though that all the cavity is filled with gas, the predicted spectro-astrometry signature of our models is a factor 2 lower than the observations by Pontoppidan et al. (2008), and the line flux is weaker than the line flux observed. Since we fit the line profile, the relative contributions to the flux by the different radii should be correct. The missing flux might be related to a CO heating mechanism or CO physics not yet included in ProDiMo.

5.3. Inner and outer disk scale height

We have set the scale height at the reference radius of 10 AU to be the same for the inner disk and the small dust grain component of the outer disk. In Fig. 8b, we display the effect of changing the scale height, the power law exponent of the surface density of the inner disk, and the gas-to-dust ratio of the outer disk in the flux of the [O i] 63 μm line.

The [O i] 63 μm line flux changes little with a change in the surface density of the inner disk, and is slightly sensitive to changes in the scale height of the inner and outer disk. In general, higher H/R values produce larger [O i] 63 μm fluxes. However, the mild sensitivity in the H/R of the outer disk is because most of the gas mass in the outer disk is in its midplane layer that has lower H/R and large grains. Changes in the scale height of the small dust extended outer disk’s atmosphere change the [O i] 63 μm line flux little because its mass is low with respect to that of the midplane.

The scale height of the disk is mainly constrained by the SED fit. In Fig. 8a, we show three insets displaying the SEDs for Models 5 with H/R = 0.7, 0.11, and 0.15. The SED is compatible with Models 5 with an H/R between 0.9 and 0.13. Disks with H/R lower than 0.9 underestimate the flux at λ> 10 μm region, and disks with H/R over than 0.13 overestimate the flux in the 10−100 μm region.

5.4. Outer disk’s gas-to-dust ratio

The fit to the SED provides a relatively robust estimate of the mass of mm-size grains in the outer disk. This value changed little in all the models tested and is on the order of 2× 10^-4M_⊙. The gas mass, hence the gas-to-dust ratio, in the outer disk is constrained by the simultaneous fit to the [O i] 63 μm line flux and the CO P(10) line profile.

Figure 8b displays the flux of the [O i] 63 μm as a function of H/R for gas-to-dust ratios of the outer disk ranging from 1 to 100. We can see that the flux of the [O i] 63 μm line is very sensitive to the value of the gas-to-dust ratio in the outer disk. We find that the [O i] 63 μm line flux tends to be described by models with gas-to-dust ratios in the outer disk that are much lower than 100. Gas-to-dust ratios between 25 and 4 provide the best fit to line flux for the scale heights that are compatible with the SED.

In all the models, the [O i] 145 μm line flux is below the Herschel flux upper limits. In all the models, CO sub-mm lines are so optically thick (τ> 100 see Fig. 6) that decreasing the gas mass does not affect the line fluxes, so these cannot be used to trace the gas mass.

The amount of gas in the inner rim of the outer disk influences the shape of the CO P(10) profile, making it more centrally peaked when more gas is present in the outer disk (see Fig. 8c). The best fit to the CO P(10) line is given by gas-to-dust ratios in the outer disk below 10.

The exact value of the gas-to-dust ratio in the outer disk is model dependent. However, since most of the models that simultaneously describe the [O i] 63 μm line flux, the SED, and the CO P(10) line profile require gas-to-dust ratios for the outer disk that are smaller than 10, we believe that a gas-to-dust ratio much lower than 100 in the outer disk is a robust result.

Finally, the simultaneous modeling of the CO P(10) line profile and [O i] 63 μm line flux favors models in which the gas surface density contrast between the inner and the outer disk at the inner rim of the outer disk is less than a factor 20.

5.5. Inner disk’s gas mass, dust mass, and gas-to-dust ratio

For a given dust composition and size distribution, the SMA 870 μm photometry upper limit of 10.5 mJy in a 34 × 70 AU beam centered on the star (Andrews et al. 2011) sets the maximum amount of dust that can be located inside 30 AU almost independently of the surface density power law exponent and the scale height of the inner disk. With the currently available data, the dust size distribution cannot be constrained.

Assuming a_min = 0.1 μm and a_max = 1000 μm for the inner disk (similar to that of the outer disk, except for a_min that is slightly larger to avoid the presence of a strong 10 μm silicate feature), and an inner disk extending from 0.2 to 30 AU, a dust mass in the inner disk of 1 × 10^-7 M_⊙ astronomical silicates reproduces the dip in the SED at 15 μm and generates a flux of 1.6 mJy at 870 μm inside 30 AU. The disk model convolved with a 34 × 70 AU beam produces 6.7 mJy of continuum flux centered on the star, which is a flux below the 3σ upper limit of the SMA observations. Using this dust mass, we found that gas-to-dust ratios between 100 and 200 can describe the CO P(10) line profile for a surface density power law exponent of 0.2. Lower gas masses in the inner disk (i.e., lower gas-to-dust ratios) generate CO P(10) profiles that are too narrow, and higher gas masses (i.e., higher gas-to-dust ratios) produce line profiles with high-velocity wings that are too broad to be compatible with the observations. The best fit was achieved with a gas-to-dust ratio of 150, hence, a gas mass inside the cavity of 1.5 × 10^-5 M_⊙.

There is a degeneracy between the gas mass of the inner disk and the exponent of the power law of the surface density. By increasing the power-law exponent (i.e., more positive), higher gas masses in the inner disk are allowed while reproducing the CO ro-vibrational line profile. However, there is an upper bound to the gas mass in the inner disk, because the surface density of the inner disk should be equal to or lower than the surface density of the outer disk at 30 AU and the flux of the [O i] line at 63 μm should be reproduced. We found that up to 10^-4 M_⊙ of gas is possible in the inner disk.

An alternative way of exploring the gas-to-dust ratio in the inner disk is to leave the gas mass and surface density constant and to decrease the dust mass in the disk at 0.2 <R< 30 AU. In Fig. 9, we display the expected CO P(10) line profiles for Model 5, keeping the inner disk’s gas mass constant (1.5 × 10^-5 M_⊙) and decreasing the amount of dust down to a factor 100 (i.e., gas-to-dust ratios from 150 up to 15 000). We found that if the gas surface density inside the cavity is constant, a lower (astronomical silicates) dust mass does not significantly affect the CO 4.7 μm profile and line flux, or the flux of the [O i] 63 μm line. A lower silicate dust mass inside the cavity leads to weaker emission in the 5−10 μm band. However, since the SED near-IR emission is dominated by the carbon grains in the innermost disk, the SED fit is still compatible with the observations. The lower bound to the dust mass inside the cavity is not constrained by current observations. To constrain the dust mass inside the cavity, observations at high spatial resolution are required, for example with ALMA. In Table 6 we predict the continuum fluxes at 430 and 870 μm inside 30 AU for different silicate dust sizes and masses in Model 5.

5.6. Outer radius of the inner disk

In the family of models 5, the inner disk has gas up to 30 AU. We have tested models with outer radii of the inner disk ranging from 15 to 30 AU. Solutions to the CO ro-vibrational line profile can be found with a similar gas mass to Model 5, by allowing for a more positive exponent of the power law of the surface density distribution. For example, for R_out = 15 AU, q increases to +0.8. Models with an inner disk outer radius smaller than 30 AU produce slightly better fits to the top of the CO ro-vibrational line; however, they produce lower spectro-astrometry signatures and generate a surface density in the inner disk’s outer radius that is slightly higher than the surface density at the outer disk’s inner radius. Current data do not allow us to set the exact value of the outer radius of the inner disk. We favor a model with a cavity replenished with gas up to 30 AU because it generates a larger spectro-astrometry signature (although even with R_out = 30 AU, the spectro-astrometry signature is underpredicted) and because surface densities with power-law exponents larger than 0.5 might result in an unstable disk. In summary, the presence of a gap in the gas of a few AU in the gas is consistent with the current data, but a large gap in the gas of tens of AU does not appear likely.

6. Discussion

6.1. The disk gas structure of HD 135344B and the origin of the transition disk shape

The growing observational evidence indicates that the gaps suggested by the SEDs and spatially resolved sub-mm observations are in fact structures that reflect the distribution of large grains, but not necessarily the distribution of small dust, and, especially the gas.

As noted by Pontoppidan et al. (2008) the presence of gas in the cavities disfavors photoevaporation as the possible origin of the transition disk shape in HD 135334B. A gas-to-dust ratio over 100 in the inner disk indicates that the mechanism depleting the material in the cavity is more efficient at depleting dust than the gas. This could potentially favor the grain-growth scenario. But, as noted by Birnstiel et al. (2012), grain-growth alone has difficulties explaining the lack of sub-mm emission in the gap.

Positive and flat power law surface density profiles (q = 0 and 1) for the inner disk of transition disks have previously been suggested by Dong et al. (2012) in the framework of their general disk model to explain the H-band scattered light present inside the sub-mm dust cavities. Our results for the modeling of the CO ro-vibrational lines in HD 135344B provide an independent argument suggesting that the surface density profile inside the cavity can increase with radius in transition disks.

What implications does a positive surface density profile have? In most protoplanetary disks the surface density distribution can be described by a power law with a negative exponent, typically −1.0 (i.e., a surface density that decreases with radius). A surface density profile with a positive power law exponent (or at least flat), indicates that the disk’s gas structure has significantly changed in the inner disk. It is not entirely clear what mechanism will generate a surface density increasing with radius.

One interesting possibility is that we are indeed observing the effects of a Jovian planet inside the inner cavity. Although the exact results of simulations depend on the inner boundary conditions set, for example the accretion rate onto the star, and the evolution time of the simulation, models studying the interaction of a Jovian planet and the disk show that a disk with an initial negative power-law surface density profile exponent could evolve into a disk that has a flat or a positive power-law surface density profile exponent at R<R_planet (see for example Fig. 2 in Lubow & D’Angelo 2006; Fig. 1 in Varnière et al. 2006; or Fig. 5 in Tatulli et al. 2011). Furthermore, a single Jovian planet is expected to open a gap typically of a few AU width (see review by Kley & Nelson 2012). Such a small gap is compatible with current data and models of HD 135344B. Another predicted effect of the presence of a planet in the disk is a higher surface density for the gas in the outer disk with respect to the surface density for the gas in inner disk. We retrieve such a gas surface density structure in our models, although, we find a surface density difference that is smaller than found by the migration calculations (e.g., Lubow & D’Angelo 2006).

Taking all this together, the characteristics found in HD 135344B evoke the effects of a migrating Jovian planet present in the disk. However, at the present time there is no firm evidence for such a companion. Our models disfavor the presence of gaps of tens of AU in the gas inside the cavity of HD 135344B, therefore, disfavor the hypothesis that a multi-Jovian planet system (e.g., Zhu et al. 2011; Dodson-Robinson & Salyk 2011) is what is responsible for the transitional disk shape in HD 135344B. This conclusion could be extended to other accreting pre-transition disks with CO ro-vibrational emission extending tens of AU or with CO sub-mm emission inside the cavity.

Finally, the result that the total the amount of gas in the disk of HD 135344B (a few 10^-3 M_⊙) is much lower than expected for a gas-to-dust ratio 100 indicates that HD 135344B is an evolved protoplanetary disk that has already lost a large portion of its gas mass. A low gas-to-dust has been reported for another transition disk, IRS 48 (Bruderer et al. 2013). Further detailed studies with a larger sample of objects are required to test whether a lower gas mass is a common characteristic of the transition disk population.

6.2. Sublimation of dust grains in the inner disk, carbon grains, and optically thick “inner walls” in pre-transition disks

One crucial ingredient for the simultaneous fitting of the SED and the CO 4.7 μm emission that extends tens of AU is the hypothesis of an innermost disk (R< 0.25 AU) enriched with carbonaceous grains. This feature is essential for making it possible for the CO warm gas at tens of AU to contribute to the CO P(10) line profile while having an inner disk that reproduces the near-IR SED.

The VLTI/PIONIER data clearly shows that the H-band excess is located inside 0.16 to 0.20 AU, the silicate sublimation radius of HD 135344B (the exact value depends on the mass of dust and n_H). This provides a plausible physical explanation for the carbon enrichment in the innermost disk: carbon grains are present at 0.08 <R< 0.2 AU because silicate grains have sublimated. The mass of carbonaceous grains required to fit the SED and the near-IR visibilities is 2.5 × 10^-12 M_⊙, which is a small amount of dust when compared to the total dust mass of the disk that is around 10^-4 M_⊙ and the dust mass at R< 30 AU that is 10^-9 to 10^-7 M_⊙. In Model 5, the carbon abundance relative to n_H required in carbon grains is around two times the solar abundance.

We assumed that the size distribution of the carbonaceous grains in the innermost disk is 0.1 to 10 μm. Future detailed work on the process of dust evaporation in the inner disk of HD 135344B (for example, as performed in the inner rim of Herbig Ae stars by Kama et al. (2009), or in the warm debris disk of Formalhaut by Lebreton et al. (2013) including the effects of gas) would be required to constrain with a physical basis the dust size distribution of carbon and silicates as a function of the radius in the inner 1 AU of HD 135344B.

The presence of carbonaceous grains in inner disks to fit near-IR visibilities has been suggested in a few circumstellar disks, such as Absil et al. (2006) and Lebreton et al. (2013) in the cases of the warm debris disks of Vega and Formalhaut. Kraus et al. (2013) analysed the (pre-)transitional disk V1247 Ori, a source that is relatively similar to HD 135344B. To fit V1247 Ori near-IR and mid-IR interferometry data, Kraus et al. suggest a disk structure where the innermost disk is composed of a mixture of carbon (50%) and silicate (50%) grains at 0.19 <R< 0.34 AU, and a mass of 10^-7 M_⊙ carbon dust grains inside 0.3 <R< 46 AU. The gap composed of carbon was suggested to avoid the presence of the 10 μm silicate feature and also to avoid a strong excess at λ> 8 μm. Furthermore, Kraus et al. suggest the presence of some inhomogeneities in the “gap” material. At the time of writing, there have been no observations of CO ro-vibrational emission for V1247 Ori. In our model of HD 135344B, we ruled out any large percentage of carbon grains (>25%) at 0.2 <R< 30 AU on the basis of the CO ro-vibrational profile. A gap with a large amount of carbon grains produces a narrow single-peaked CO ro-vibrational line profile because the CO emission is dominated by the inner rim of the outer disk (see Model 2).

The carbonaceous grains innermost disk inside the silicate sublimation radius proposed here for HD 135344B provides an alternative disk structure to the “optically thick-wall” suggested for pre-transition disks (e.g., Espaillat et al. 2007, 2010; Brown et al. 2007). We have shown here that the “optically thick wall disk structure” is not compatible with detecting CO ro-vibrational emission extending tens of AU. Therefore, the carbon-rich (or refractory grains rich) inner disk structure suggested for HD 135344B might be applicable to other transition disks displaying near-IR excess and CO 4.7 μm emission extending several AU (e.g., SR 21, Pontoppidan et al. 2008). Further near-IR interferometry observations of (pre-)transition disks are required to establish what fraction of them have dust material inside the silicates sublimation radius. Moreover, since CO ro-vibrational emission has been detected in several transition and pre-transition disks (e.g., Pontoppidan et al. 2008; Salyk et al. 2009), it would be of great interest to extend the analysis presented here to them.

The models described here have used amorphous carbon grains. We tested graphite grains using a mixture of 50% parallel and 50% perpendicular opacities. The fit to the CO P(10) line profile did not change significantly. The SED displayed a slightly weaker emission at 3−10 μm.

A potential challenge for the presence of carbonaceous grains in the innermost disk (0.08 <R< 0.2 AU) is an oxygen-rich atmosphere. Carbon grains would react with the oxygen and be destroyed (see for example the models of Gail 2001; Lee et al. 2010). The presence of carbon grains in the innermost disk depends on the balance between the destruction and the replenishment timescales.

To estimate the destruction timescale by oxygen atoms (i.e., chemical sputtering), we take one carbon grain of density ρ, initial size a₀, and radius a. Atomic oxygen reacts with individual carbon atoms at the surface of the carbon grain, producing a change in the carbon grain’s size with time (Evans 1994): $\frac{\mathrm{\Delta }\mathit{a}}{\mathrm{\Delta }\mathit{t}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{0}}}{\mathit{t}}\mathrm{\simeq }\mathit{ȧ}\mathrm{=}\mathrm{-}\frac{{\mathit{n}}_{\mathrm{O}}{\mathit{m}}_{\mathrm{C}}\mathit{Y}}{\mathit{\rho }}{\left(\frac{\mathit{k}{\mathit{T}}_{\mathrm{gas}}}{\mathrm{2}\mathit{\pi }{\mathit{m}}_{\mathrm{O}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(2)where t is the time for the grain to be destroyed, n_O the number density of oxygen (n_O = n_Hχ_O and χ_O = 3 × 10^-4), Y the yield of this Eley-Rideal type reaction, k the Boltzmann constant, and m_O the mass of the oxygen atom (m_O = 16 × m_H). Here we assume that the destruction rate is constant. Solving for t we obtain $\mathit{t}\mathrm{~}{\mathit{a}}_{\mathrm{0}}\frac{\mathit{\rho }}{{\mathit{n}}_{\mathrm{O}}{\mathit{m}}_{\mathrm{C}}\mathit{Y}}{\left(\frac{\mathrm{2}\mathit{\pi }{\mathit{m}}_{\mathrm{O}}}{\mathit{k}{\mathit{T}}_{\mathrm{gas}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(3)If the grain’s mass density ρ is 2.2 g cm^-3, then the life time of a carbon grain of radius a₀ is $\mathit{t}\mathrm{~}\mathrm{0.28}\mathrm{\times }\left(\frac{\mathrm{0.1}}{\mathit{Y}}\right)\left(\frac{{\mathit{a}}_{\mathrm{0}}}{\mathrm{0.01}\mathit{\mu }\mathrm{m}}\right){\left(\frac{\mathrm{2000}\mathrm{K}}{{\mathit{T}}_{\mathrm{gas}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{{\mathrm{10}}^{\mathrm{10}}{\mathrm{cm}}^{-3}}{{\mathit{n}}_{\mathrm{H}}}\right)\mathrm{yr}\mathit{.}$(4)The value of the yield is not well constrained. Assuming Y = 0.1 (Vierbaum & Roth 2002, but it can be lower), n_H = 2 × 10¹⁰ cm ^-3 (midplane value at 0.1 AU in Model 5, see Fig. A.4), this simple calculation shows that carbon grains with sizes 0.01, 0.1, and 1 μm survive at 2000 K time scales of around a month, a year, and ten years, respectively.

To calculate the carbon replenishment time scale at R< 0.2 AU, three input parameters are needed: the gas mass accretion rate, the gas-to-dust ratio, and the fraction in mass of dust grains in the form of carbonaceous grains at R> 0.2 AU. Sitko et al. (2012) suggest a gas accretion rate of 2 × 10^-8 M_⊙/yr. Our model indicates that the gas-to-dust ratio at R< 30 AU is 150 (it can be higher). The mass fraction of carbon grains with respect to silicates at 0.2 <R< 30 AU can be constrained observationally by the maximum mass of carbon that can be added and still keep the fit to the SED and the CO ro-vibrational line. We have found that up to 3% of the dust mass (10^-9 M_⊙) at 0.2 <R< 30 AU can be in the form of carbon grains and still fit the SED and CO lines. Taking the gas accretion rate, the gas-to-dust ratio, and the fraction in mass of carbon together, we obtain a supply rate of carbon grains to the innermost disk of 4 × 10^-12 M_⊙/yr. The life time of carbon grains larger than 0.1 μm is years, thus, the disk can have around 10^-12 M_⊙ transient carbon grains at R< 0.2 AU. A carbon mass similar to that is required to fit the near-IR SED in our models.

In summary, taking the supply rate via dust mass accretion and the destruction timescales in the innermost disk together, we find that a carbon innermost (R< 0.2 AU) disk could indeed be plausible in HD 135344B.

In a recent work, Lee et al. (2010) have studied the survival of carbon grains in a T Tauri disk. They find that carbon grains can be destroyed by oxygen atoms in the warm (T ~ 500 K) upper layers of the disk on time scales ranging from tens to thousands of years depending on the carbon grain size, porosity, and distance of the star. The potential destruction of carbon grains at R> 0.2 AU could present a challenge to the interpretation presented here. Lee et al. models are not dynamical models in the sense that the dust components are not followed with the gas in a simulation. Thus, it is unknown what fraction of the total amount of carbon grains will effectively reach the upper layers of the disk and be destroyed, what fraction of the carbon grains will survive in the midplane, and with what rate carbon grains at a few AU disk are replenished by the outer disk in an accreting disk. Future, more detailed models that include the dynamics of the gas, the interaction between the gas and the dust, the photochemistry of the disk, and carbon destruction would be of great help in further testing of the scenario proposed here.

From the point of view of interpreting the observations, other refractory grains in the innermost disk, such as titanium, calcium, or aluminum-bearing minerals (e.g., titanium or aluminum oxide), could be responsible for the near-IR excess observed. Since the current data do not constrain the inner disk composition to this degree of refinement (i.e., multiple refractory dust species), we did not attempt to fit our data with mixtures of carbon and/or other refractory components. We retained the amorphous carbon solution since it is the simplest refractory grain that can be assumed.

The conclusion that the innermost disk of HD 135344B could be carbon-rich potentially presents a problem in a wider cosmochemical context, because the Earth and other terrestrial planets (as well as some exoplanetary systems, e.g, Xu et al. 2013) display a carbon to silicon deficit in comparison with meteorites, the Sun, and the Interstellar Medium. Perhaps HD 135344B is an unique object, or the processes responsible for the carbon depletion on Earth are particular to the Solar System, or perhaps there is much wider variety of the possible chemical compositions that terrestrial planets can have in different planetary systems. We stress that the carbon-rich innermost disk of HD 135344B is a plausible solution that allows us to simultaneously reproduce all available gas and dust observables of the inner disk, most notably the near IR visibilities indicating continuum emission inside the silicate sublimation radius and the CO ro-vibrational emission extending tens of AU. But other refractory grains could provide similar solutions. Amorphous carbon grains provide the simplest solution that fits the data, but they are not the only possible solution.

6.3. Double- and single-peaked CO-rovibrational line profiles in transition disks

In a recent study of CO-rovibrational emission in young stars, Bast et al. (2011) have observed three types of line profiles: (1) double-peaked; (2) narrow-single peaked; and (3) single-peaked with a broad base. Bast et al. argue that the profiles of type 3 stem from a combination of emission from the inner part of the disk (R< a few AU), and emission from a slow moving disk wind. In their analysis, Bast et al. suggest the line profile parameter P₁₀ (the full width line profile at 10% of its height divided by the full width of the line profile at 90%) to distinguish between the different types of line profiles. In the case of HD 135344B, the value of P₁₀, obtained from an average of the lower J-transition lines up to P(14), is 6.9. Thereby locating HD 135344B just above the value of 6, the maximum line profile parameter for a Keplerian disk model with a power-law temperature profile.

In our models of HD 135344B, we have found that the CO ro-vibrational line profile is the result of the combination of two contributions: (1) the emission from the gas inside the dust cavity and (2) the emission from the gas located in the inner rim of the outer disk. The relative contribution of both components, hence, the CO ro-vibrational line profile, depends on the properties of the dust in the inner disk and the gas mass in the inner and outer disk.

In Models 1 and 2, we have shown the effect of the inner disk dust composition on the CO ro-vibrational line profile. If the optical depth of the continuum at 4.7 μm inside 30 AU is small (τ< 0.1, e.g., Model 2), then the CO ro-vibrational emission will be dominated by the contribution of the inner rim of the outer disk, and a narrow single-peaked line profile would be produced.

A single-peaked line profile can also be obtained by increasing the gas mass in the outer disk, diminishing the gas mass in the inner disk, or increasing the scale height of the outer disk with respect to the inner disk. In the case of Model 5, if the gas mass in the outer disk is increased such that the gas-to-dust ratio is more than 10, then the line profile becomes more and more single-peaked (see bottom panel of Fig. 8).

In summary, our radiative transfer calculations of the CO 4.7 μm emission show that in the case of transition disks, the presence of single-peaked CO ro-vibrational line profiles can be explained by Keplerian disk emission without recourse to a disk wind component. This only applies to sources in which the emission line center is not blueshifted with respect to the stellar velocity.

7. Summary and conclusion

We conducted a modeling project aimed at constraining the gas mass and the gas and dust disk structure of the transition disk HD 135344B from multi-instrument and multiwavelength observations of gas and dust. We found that the previously suggested inner disk structure (Brown et al. 2007), namely a narrow dust inner disk of from 0.18 to 0.45 AU followed by a large 45 AU dust gap replenished with gas, fails to reproduce the CO ro-vibrational emission observed as the line profile produced from this disk model is a broad double peak.

We have found a disk model that is able to reproduce current observational constraints. This disk is composed of three zones:

• 1.

  A first zone between 0.08 and 0.2 AU composed of small carbonaceous grains (and gas) with a total dust mass of a few 10^-12 M_⊙ (a few solar abundances of carbon). The presence of these inner carbonaceous grains provides

  • a)

    a fit to the near-IR H -band visibilities and closure phases,

  • b)

    a fit of the near-IR SED, while allowing the warm CO at several AU to emit and contribute to the 4.7 μ m line profile,

  • c)

    an agreement with the higher temperatures ( T> 1500 K) expected in this zone.

• 2.

  A second zone extending from 0.2 to 30 AU (i.e., the dust cavity) replenished with gas (10^-5–10^-4 M_⊙) with a surface density increasing as a function of the radius and dust mass of astronomical silicates of maximum 10^-7 M_⊙. An increasing surface density profile with radius is required to fit the shape of the CO ro-vibrational emission lines. The fit to the SED constrained the scale height between 0.09 and 0.13 at 10 AU with a flaring exponent 1.12. The gas-to-dust ratio in this zone is larger than 100, however, the exact value is not well constrained. We found models up to gas-to-dust ratios of 15 000 that are consistent with the observations either by decreasing the silicate dust mass by a factor 100, or by increasing both the gas mass by a factor of a few and the power-law exponent of the surface density distribution. The upper bound to the gas mass at R< 30 AU is given by the flux of the [O i] 63 μm, combined with the requirement that the surface density of the inner disk should be equal or lower than the surface density of the outer disk at 30 AU. The dust surface density at R< 30 AU is lower than the one expected from extrapolating the dust surface density from the outer disk. This zone can contain up to a 3% in mass (~10^-9 M_⊙) in carbon grains and keep the fit to the SED.

• 3.

  A third zone from 30 AU to 200 AU (the outer disk) with astronomical silicate grains, a dust mass of 2 × 10^-4 M_⊙, a gas mass 10^-4−10^-3 M_⊙ (gas-to-dust ratio <10), surface density exponent of −1.0, and flaring of 1.0. In this zone large (0.1 <a< 1000 μm) and small (0.1 <a< 10 μm) dust grains have different radial and vertical spatial distributions. Large grains are located at 45 <R< 200 AU in a disk with scale height 0.05 and 75% of the gas and dust mass. Small grains are located at 30 <R< 200 AU in a disk with higher scale height (0.09 to 0.13, i.e., the same H/R as zone 2) and 25% of the gas and dust mass. The vertical structure of the outer disk echoes the expected effects of dust growth and settling. A gas-to-dust ratio much lower than 100 in the outer disk is required to fit the [O i] 63 μm line flux and to reproduce the CO ro-vibrational line profile simultaneously.

The models suggest that the best fit to the gas observations in HD 135344B is provided by a disk in which the gas surface density and the scale height have no large discontinuities at 30 and 45 AU. In other words, there is no large gap in the gas distribution of HD 135344B. The cavity observed in the near-IR and sub-mm is replenished by gas and (some) dust. The presence of a small gap of a few AU in the gas is consistent with current data, a large gap in the gas of tens of AU does not appear likely.

The global gas-to-dust ratio, i.e., integrated over the full disk, is much lower than 100. This provides further evidence that HD 135344B is an evolved protoplanetary disk that has already lost a large amount of its gas mass. The disk structure proposed for HD 135344B could be applied to other pre-transition disks with CO ro-vibrational emission extending to several AU.

The increasing surface density profile of the gas in the inner disk, the difference in the radial distribution of large and small grains in the outer disk, and a small gap in the gas between the inner and outer disk are compatible with the dynamical interaction of a single Jovian planet and the disk. However, we should be cautious because other mechanisms could be responsible for the gas and dust distribution observed in HD 135344B.

The current ensemble of observations show that the transitional disk features (i.e., gap) observed in the continuum reflect only changes on the distribution of large and small dust particles in the disk. The gas has a spatial distribution that can differ from that of the dust, particularly that of large grains observed in the sub-mm. Furthermore, we find that the gas distribution and mass in transition disks can be very different from primordial disks.

Our modeling results predict that gas and (some) dust emission should be detected inside the cavity with high-sensitivity and high-spatial resolution sub-mm observations. ALMA measurements of continuum emission inside 30 AU would be a great help in determining the dust size distribution and dust mass inside the cavity, quantities that at present are free parameters in the model. ALMA gas observations at high spatial resolution, such as the CO 6−5 line, used in conjunction with other gas tracers will enable the gas mass and surface density inside the cavity and in the outer disk to be better constrained.

Future detailed dynamical disk models that include the dynamics of the gas, the interaction between the gas and the dust, the photochemistry of the disk, and carbon destruction would be a big help for establishing whether the carbon innermost disk suggested by our modeling is viable or whether including refractory grains different from carbon is necessary.

In this paper we showed that the simultaneous modeling of gas and dust observations is required to address the problem of the dust and gas structure of protoplanetary disks.

Online material

Appendix A

A.1. Photometry

In this appendix, we summarize the different sources of the photometry used to construct the SED (see Table A.1).

A.2. Stellar parameters

The stellar properties of HD 135344B have been studied by previous authors (e.g., Dunkin et al. 1997; Müller et al. 2011; Andrews et al. 2011). While authors agree on a spectral type F4, an effective temperature 6750 ± 150 K (6660 K, Dunkin et al. 1997; 6590 K, Andrews et al. 2011; 6810 ± 80 K, Müller et al. 2011), and a mass of 1.6 to 1.7 M_⊙, there is a discrepancy in estimations of the radius of the star. Müller et al. (2011) suggest a radii of 1.4 ± 0.25 R_⊙, while Andrews et al. (2011) suggest a radii of 2.15 R_⊙.

To have an independent assessment of the star properties, we downloaded reduced archival high-resolution (R = 75 000) ESO/VLT-UVES spectra of HD 135344B in the 4800 to 5500 Å and 5830−6800 Å ranges and used the interactive spectra visualization tool described in Carmona et al. (2010) to compare the UVES spectrum to BLUERED (Bertone et al. 2008) high-resolution synthetic spectra, in regions not contaminated by emission lines or telluric absorption. We found that the optical spectra of HD 135344B is compatible with spectral templates with T_eff ranging from 6500 to 7000 K for logg ranging between 4.0 and 5.0 (see Fig. A.1). Naturally, logg = 5.0 is not realistic as logg = 4.33 for an F-type star on the main sequence. The minimum in the χ² statistic suggests a T_eff around 6750 K.

To have a second constraint on the spectral type and A_V, we employed the HST/COS⁶ spectrum of HD 135344B. In Fig. A.2, we display coadded COS spectrum of HD 135344B, smoothed by a running nine-point boxcar filter, and dereddened by a Cardelli et al. (1989)R = 3.1 extinction curve with E(B − V) = 0.129 (Av = 0.4). We overplot F3V (HR 9028, IUE SWP 14002, red) and F4V (HD 27901, SWP 45935, blue) spectral type comparison stars scaled by the difference in V magnitudes between them and HD 135344B. The spectral template F4V provides the best fit to the observed COS spectrum. The F4V spectrum corresponds to a T_eff of 6620 K using Schmidt-Kaler (1982).

To obtain an estimate of the luminosity and to constrain the value of the radius, we used a distance of 140 pc (van Boekel et al. 2005), the B and V photometry (9.2 and 8.7, SIMBAD database), and the T_eff, M_V, BC values of Schmidt-Kaler (1982) for stars of luminosity class V (i.e., smallest possible radii). We found that for a T_eff equal to 6440 K (F5) and 6890 K (F2), R = 2.09 and 2.07 R_⊙ respectively. Interpolating for T_eff 6620 K, we obtained 2.08 R_⊙. These values are closer to the radius estimate of 2.15 R_⊙ suggested by Andrews et al. (2011) than the 1.4 R_⊙ claimed by Müller et al. (2011).

For our models, we used a star with T_eff = 6620 K (F4V), a stellar radius of 2.1 R_⊙, a mass of 1.65 M_⊙, and A_v = 0.4.

A.3. CO 4.7 μ m emission

A.4. Model 5: Herschel CO and water lines

A.5. Details of Model 5

In this figure, we show additional plots of the structure and emitting regions of several gas tracers in Model 5.

---

Acknowledgments

A. Carmona and C. Pinte acknowledge funding from the European Commission’s 7th Framework Program (EC-FP7) (contract PERG06-GA-2009-256513) and from the Agence Nationale pour la Recherche (ANR) of France under contract ANR-2010-JCJC-0504-01. Calculations were performed at the Service Commun de Calcul Intensif de l’Observatoire de Grenoble (SCCI) on the Fostino supercomputer funded by the ANR (contracts ANR-07-BLAN-0221, ANR-2010-JCJC-0504-01 and ANR-2010-JCJC-0501-01) and the EC-FP7 contract PERG06-GA-2009-256513. The research leading to these results has received funding from the EC-FP7 under grant agreement no 284405. F.M. acknowledges support from the Millennium Science Initiative (Chilean Ministry of Economy), through grant Nucleus P10-022-F. PIONIER is funded by the Université Joseph Fourier (UJF), the Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), the Agence Nationale pour la Recherche (ANR-06-BLAN-0421 and ANR-10-BLAN-0505), and the Institut National des Science de l’Univers (INSU PNP and PNPS). The integrated optics beam combiner is the result of a collaboration between IPAG and CEA-LETI based on CNES R&T funding. Part of this work was funded by the Agence Nationale pour la Recherche of France through the Chaire d’Excellence grant ANR (CHEX2011 SEED). This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France. This work is based (in part) on archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. We thank Alex Brown for providing the Chandra X-ray luminosity.

References

All Tables

All Figures