© ESO, 2015

1. Introduction and motivation

It is commonly assumed that planets form in circumstellar disks around young pre-main-sequence (PMS) stars (e.g., Mordasini et al. 2010). During this formation process planets perturb the disk because of their gravitational potential (Goldreich & Tremaine 1980). The results of this interaction are characteristic structures such as spiral arms (Ward 1997; Tanaka et al. 2002) and gaps (Papaloizou & Lin 1984). Therefore, observations of these large-scale structures provide indirect constraints on the existence of protoplanets and on the relevant physical processes involved in the protoplanet-disk interaction.

Wolf & D’Angelo (2005) and Ruge et al. (2013) showed the possibility of tracing gaps in continuum re-emission light for various disk configurations and ALMA configurations. In addition, Ruge et al. (2014) investigated the feasibility of tracing gaps in continuum re-emission maps and scattered light images simultaneously. In contrast to this, only a few studies exist investigating the feasibility of tracing planet-induced structures through molecular line observations. Semenov et al. (2008) showed the feasibility of drawing conclusions about disk properties such as large-scale temperature gradients and chemical stratification from simulated HCO⁺ (4−3) observations. Dutrey et al. (2008) used spatially unresolved but spectroscopically resolved line observations to trace inner cavities due to the radial dependence of the (Keplerian) velocity field. Cleeves et al. (2011) demonstrated the feasibility of detecting a huge inner cavity (R ~ 45 AU) using molecular lines (CO, H₂CO₂). Using ALMA, Casassus et al. (2013) were able to image multiple molecular gas species within the dust free inner hole of HD 142527 and to utilize the varying optical depths of the different species to detect a flow of gas crossing the gap opened by a planet.

However, the main component of protoplanetary disks, cold molecular hydrogen (H₂), is very difficult to observe directly in these disks because of the missing electrical dipole moment. Thus, other molecules with lower excitation temperatures must be observed to derive the dynamics, temperature, density, and chemical structure in these disks. The most frequently observed molecule is CO, mainly because of its high abundance and low excitation energy. Observations of various isotopologues of this molecule provided the first constraints on the physical properties such as the (vertical) gas temperature structure of young disks (e.g., Piétu et al. 2007). Through comprehensive surveys of protoplanetary disks around solar-type stars other molecules such as HCO⁺, H₂CO, CS, CN, and HCN were detected (Kastner et al. 1997; Öberg et al. 2010; Thi et al. 2004). Recently, Chapillon et al. (2012a) provided evidence of HC₃N in the disks around GO Tau and MWC 480 using the IRAM 30 m telescope. However, because of insufficient spatial resolution and the lack of sensitivity these observations were restricted to the simplest and most abundant molecules. The Atacama Large (Sub)-Millimeter Interferometer (ALMA) will potentially allow the direct study of the molecular composition and distribution in protoplanetary disks in the submillimeter wavelength regime both at global scale as well as in the potential planet-forming region. The first ALMA observations of HD 163296 show that CO radiation originates from a thin layer with an opening angle of about 15° with respect to the disk midplane (Rosenfeld et al. 2013; de Gregorio-Monsalvo et al. 2013). More recently with the unprecedented sensitivity of ALMA, cyclic C₃H₂ (which has been proposed as a useful probe of radiation penetration in disks) has been detected in the disk surrounding HD 163296 (Qi et al. 2013).

The motivation to search for gaps in protoplanetary disks using molecular line observations is given not only by the aim to indirectly detect planets, but even more to derive complementary constraints on the disk physics from the planet-disk interaction. For example, the resulting gap structure, given by its shape and depth in the gas and dust distribution, depends not only on the mass of the planet and grain size, but also on the interaction of the planet with the gas and dust phase (Paardekooper & Mellema 2004). As a result, larger dust grains (≳150 μm) are expected to decouple from the gas, pile up at the edges of the gap and significantly alter the gas-to-dust mass ratio (Paardekooper & Mellema 2006). Moreover, we emphasize that a gap in the density distribution of a cirumstellar disk is not necessarily an irrevocable sign of an embedded planet. Flock et al. (2015) showed that MRI turbulence is sufficient to generate a gap in the gas phase without the presence of a planet. Recently, Bruderer et al. (2014) were able to detect CO emission inside the dust cavity of IRS 48. They concluded that a massive inner companion is the main driver for the clearing of the cavity and excluded photoevaporation and grain-growth only. Therefore, complementary observations of the gas and the dust phase are crucial in order to improve and verify our understanding of the predominant physics; e.g., the interaction between the gas and dust component in the planet forming regions of protoplanetary disks. As a first step into this potentially new field of studying disk physics through planet-disk interaction, we investigate the feasibility of tracing planet-induced gaps through spatially and spectroscopically resolved molecular line observations. For this purpose we developed the N-LTE line and dust continuum radiative transfer code Mol3D, which is based on the radiative transfer code MC3D (Wolf et al. 1999; Wolf 2003).

Our studies are based on spatial density and velocity distributions resulting from 3D hydrodynamic simulations of planet-disk interaction using the PLUTO code (Mignone et al. 2007). The temperature structure and corresponding molecular line maps are calculated with Mol3D. These velocity-channel maps are used as an input for the CASA (ver. 4.2) simulator (Petry & CASA Development Team 2012). Finally, we investigate the feasibility of tracing planet-induced gaps in the simulated ALMA maps. Based on these results, our goal is to prepare future ALMA observations of similar protoplanetary disks. In contrast to earlier studies, which were often focused on exemplary case studies (e.g., Cleeves et al. 2011), we explore a large parameter space covering various disk properties such as their total mass, size, and abundant molecules of typical PMS T Tauri and Herbig Ae stars.

This study is structured as follows. In Sect. 2 we briefly present the radiative transfer code Mol3D, discussing its underlying physical assumptions and test cases. In Sect. 3 we describe the study of tracing gaps using molecular lines with ALMA. Finally, in Sects. 4 and 5 we discuss our results and provide conclusions about the feasibility of tracing planet-induced gaps though molecular line observations with ALMA.

2. Molecular line radiative transfer code

In this section, we briefly describe our line and dust continuum radiative transfer code Mol3D.

Mol3D features

The N-LTE 3D parallelized line radiative transfer code Mol3D features full 3D spherical, cylindrical, and Cartesian grids with several types of boundary spacing (e.g., logarithmic, linear). Mol3D typically uses three steps to produce velocity-channel maps, spectra, dust continuum maps, and spectral energy distributions (SED). In the first step, the dust temperature is calculated based on the assumption of local thermal equilibrium. To achieve a high efficiency, we use a combination of the continuous absorption method proposed by Lucy (1999) and immediate temperature correction by Bjorkman & Wood (2001).

In the second step the level populations are calculated. As discussed in Appendix B.2.1, this is the main problem of line radiative transfer. In the case of circumstellar disks, the problem can be significantly simplified by using adequate approximation methods. Mol3D allows the local thermodynamic equilibrium (LTE), the free-escape probability (FEP), or the large velocity gradient (LVG) method to be applied. It has been shown that LTE is a reliable assumption for lower transitions of abundant molecules, but it is recommended that LVG be used, especially for high transitions of complex molecules with low abundances (Pavlyuchenkov et al. 2007).

The last step is to generate velocity spectra or velocity-channel maps using a new efficient ray-tracing algorithm. Here, we have to consider that the source function is not constant inside a single cell, which is due to the underlying velocity field.

As our code solves the radiative transfer equation including both gas and dust, we have to take the optical properties of both into account. For the intensity integration, we implemented an embedded Runge-Kutta-Fehlberg solver of order 4(5) with automatic step control. The main advantage is that with providing a relative (10^-8) and an absolute (10^-20) error limit for the intensity integration, the algorithm calculates the corresponding step-width automatically. Conservation of energy is ensured by a recursive pixel refinement algorithm based on the underlying model grid.

Mol3D is tested and benchmarked against other line and continuum radiative transfer codes. The results are presented in Appendix B along with a detailed derivation of the underlying concepts and physical assumptions. The code will soon be publicly available¹ under an open-source license.

3. Model setup

In this section we investigate the feasibility of observing planet-induced structures using molecular lines with ALMA. For this purpose we perform 3D hydrodynamic (HD) simulations, follow-up line-radiative-transport (LRT) calculations, and simulate ALMA observations using the CASA simulator.

3.1. Three-dimensional hydrodynamic simulations

In preparation for this study we started with an analytic density distribution including a radially symmetric sharp gap completely depleted of dust and gas. We find that this approach is insufficient because it is based on an unrealistically high density contrast between the gap and its surroundings, resulting in a gap detection rate of almost 100% for the planet-disk configurations considered in our study (see Sect. 3.2). Instead, justified assumptions about the ratio between the density inside and outside the gap, the gap width, and the depth of the gap are needed.

To reach this goal, we make use of full 3D hydrodynamic simulations obtained with the well-known finite-volume fluid dynamics code PLUTO (Mignone et al. 2007). The code uses the HLLC² solver for the hydrodynamics, a second-order Runge-Kutta time integration and second-order variation diminishing (TVD) scheme. The setup of the simulations was originally presented in Uribe et al. (2011) and the models used in this work are shown in Ruge et al. (2013). From their entire set of planet-disk configurations we consider two scenarios. The first model includes an embedded Jupiter-mass planet (fixed mass ratio M_planet/M_star = 0.001) on a fixed circular orbit. The second disk with the same initial parameter setting does not contain a planet, and thus serves as an undisturbed smooth reference disk for comparison.

We use a spherical grid where θ is the angle from the vertical axis and φ is the azimuthal angle. To save computation time, the vertical extent of our PLUTO simulations is chosen to cover only θ∈ [π/2−0.3, π/2 +0.3], covering four pressure scale heights above the disk midplane.

As initial conditions, it is assumed that the disk gas has a sub-Keplerian rotation profile, and that the density profile ρ₀ ∝ r^{− a} and the speed of sound c_s = c₀(rsinθ)^{− b} follow a power-law with exponents a = 3/2 and b = 0.5, respectively. The initial density distribution is given by (1)The disk is described by a locally isothermal equation of state . The ratio of the pressure scale height h to the radial coordinate of the disk is taken to be a constant such that h = H/ (rsinθ) = 0.07. A softened gravitational potential is used for the planet, (2)where ϵ is the softening parameter which is defined as a fraction of the Hill radius ϵ = l·r_p(M_p/ 3)^1/3 with l = 0.3.

The original disk configuration has an inner and outer edge of 2 AU and 9 AU, respectively. The planet is located at r_p = 5 AU. For the hydrodynamics simulations, buffer zones are included to avoid boundary effects, so the complete radial domain is from 1 − 10 AU. For the follow-up radiative transfer, we need to convert the resulting density distribution. We average every two cells in each direction to save memory. Thus, we are able to reduce the required number of grid cells to . For the subsequent radiative transfer (RT) simulations with Mol3D we add two cells (one above and one below the disk) to fill the entire θ∈ [0, π] domain. The grid parameters for both codes are summarized in Table 1.

3.2. Disk model

The hydrodynamic simulations discussed in Sect. 3.1 rely on arbitrary units for length and mass. This allows us to introduce linear scaling parameters to scale the disk mass and size to our needs for protoplanetary and/or transitional disks.

Disk size

For the inner and outer disk radius we chose the linear scaling parameter k ∈ [7,..,25]. Hence, the inner radius of the disk is in the range of R_in = k·2 AU (i.e., R_in∈ [14,...,50 AU]) and the outer radius in the range of R_out = k·9 AU (i.e., ∈ [63,...,225 AU]. As a result, our disk models mimic different evolutionary states from early protoplanetary disks to more evolved transitional disks or circumbinary disks with huge inner cavities. The disk density distribution resulting from the simulation of a Jupiter-mass planet located at 80 AU in a disk ranging from R_in = 32 AU to R_out = 144 AU (corresponding to a scaling parameter of k = 16) is presented in Fig. 1. Several specific features in the disk density distribution that result from planet-disk interaction can be seen. The most prominent one is the distinct gap between 68 and 92 AU. In the cut through the density distribution perpendicular to the disk midplane, the gap appears as a tunnel in the protoplanetary disk, as the planet has hardly any influence on the disk surface (see Fig. 1, right). Furthermore, the midplane density distribution (Fig. 1, left) shows pronounced spiral arms; these features are barely visible in the surface of the disk in this case as well.

Mass range

We consider a wide range of disk masses (gas+dust) from 2.67 × 10^-1 M_⊙ to 2.67 × 10^-5 M_⊙ which approximately covers the range of protoplanetary and transitional disk masses derived from submm/mm observations (e.g., Andrews et al. 2010, 2011).

Dust properties and stellar parameters

For the dust we use a mixture of 62.5% astronomical silicate and 37.5% graphite (Weingartner & Draine 2001). Furthermore, we assume the dust grains to be spherical following a grain size distribution given in Eq. (3) with an exponent of d = 3.5 (Dohnanyi 1969): (3)We consider dust grains with a minimum radius of a_min = 5 nm and a maximum radius of a_max = 250 nm as found in the interstellar medium (Mathis et al. 1977). The radiation source is the central star. For this study we consider a T Tauri and a Herbig Ae star with parameters given in Table 2).

Distribution of the different molecular species

For the gas-phase, we assume a gas-to-dust-mass ratio of 100:1 and perfect thermal coupling between gas and dust.

For the spatial distribution of different molecular species relative to each other only weak constraints are available from observational and numerical studies. For example, observations have shown carbon monoxide emission from the cold outer parts of protoplanetary disks where CO should be frozen out (<~20 K; Dartois et al. 2003; Goto et al. 2012). A possible explanation has been proposed by Aikawa & Nomura (2006) and Semenov et al. (2006). They argue, that the CO abundance could be significantly higher in the cold (midplane) regions of protoplanetary disks owing to vertical and radial mixing. Motivated by this idea, Hersant et al. (2009) computed the chemical evolution of circumstellar disks. They included grain surface reactions with and without turbulent mixing and CO photodesorption. However, the resulting CO column densities are still not consistent with observations. As indicated by Cleeves et al. (2011), the disk midplane in the gap region can be irradiated by the stellar radiation and X-rays. Thus, CO would not be frozen out onto grains, which would result in better gap detection rate thanks to higher contrast.

Because the debate on relative molecular abundances in circumstellar disks is still open, we decided to use fixed abundances relative to molecular hydrogen. In particular, we follow Pavlyuchenkov et al. (2007) and assume their reported abundances for all molecules considered in this study (see Table 3). These values are obtained by using a gas-grain chemical model with surface reactions (Semenov et al. 2004, 2005), which is based on the UMIST 95 database of gas-phase reactions (Millar et al. 1997). This approach keeps the model simple and hence reduces the complexity of our analysis.

Line radiative transfer setup

The level populations are calculated by using the LVG method (see Appendix B.2.1) which has proven its reliability for protoplanetary disks (Pavlyuchenkov et al. 2007). Using these results, we produce synthetic velocity-channel maps and spectra for every model configuration. We assume a pure Keplerian rotation profile. In this study we fix the turbulent velocity v_turb = 100 m/s for all of our models, following the results of earlier molecule observations (e.g., Piétu et al. 2007; Chapillon et al. 2012b). We choose a fixed bandwidth of 50 m/s and calculate 141 velocity-channel maps for each of the models with 10° inclination and 71 for every face-on orientated model in our sample. An overview of the considered parameter space of our model is compiled in Table 4.

4. Results

4.1. Ideal line maps

In this section we present ideal velocity-channel maps, analyze the appearance of planet-induced disk perturbations, and investigate the feasibility of tracing these structures with ALMA (Sect. 4.2). In the following we use the term “depth of the gap” as the flux density contrast in the velocity-channel maps caused by the gap.

Optical depth effects

We begin our discussion based on ideal ¹²C¹⁶O (3−2) velocity-channel maps of our T Tauri disk model. To start with, we assume a total disk mass of 2.67 × 10^-5 M_⊙, which is at the very low end of typical disk masses of these

objects (e.g., HD 166191, Coka Tau/4; Kennedy et al. 2014; Forrest et al. 2004; D’Alessio et al. 2005). The inner rim of the chosen disk has a radius of 32 AU while the outer radius amounts to 144 AU. The planet is located at ~80 AU and the planet-induced gap has a width of ~32 AU. The disk is seen face-on. In this case, the line shows a symmetric gap throughout the different velocity channels, because the line width is given only by thermal and microturbulent velocities (see Fig. 2; see also Appendix B.2). We find that the contrast between the gap and its edges is not constant throughout the individual velocity channels, but varies with the optical thickness of the chosen line. In the optically thin case, the gap depth increases as it approaches the line center owing to the increasing line flux. In our case, the ¹²C¹⁶O (3−2) transition is optically tick at its center even for this lowest disk mass of our sample. Thus, the τ = 1 surface seen in the line wings belongs to parts of the disk closer to the disk midplane than the disk layers traced at the line center (i.e., at 0 m/s). Therefore, for this specific configuration we observe the gap best at velocities ±300 m/s, while the main contribution at the line center comes from the upper layers of the disk. As these upper layers are less disturbed by the planet-disk interaction, the gap is less pronounced in the center of an optically thin line, i.e., the nearly undisturbed upper disk layers mask the gap. In analogy, further signatures of the planet-disk interaction like the pronounced spiral arm in the outer disk can be identified best in the line wings (here, in the ±400 m/s channel maps).

If we increase the optical depth by increasing the disk mass by two orders of magnitude (i.e., 2.67 × 10^-3 M_⊙), which represents a typical disk mass of protoplanetary disks, we find that the gap vanishes nearly completely in the velocity-channel maps (Fig. 3). As only the envelope above the gap can be observed, the disk appears to be undisturbed without any clear signs of planet-disk interaction. This is the main problem when observing ¹²C¹⁶O. For disk masses corresponding to typical protoplanetary disks masses we find that the ¹²C¹⁶O molecule is not an adequate tracer for gaps. In this case, other less abundant molecules/isotopologues that allow deeper disk layers to be observed are better suited. We illustrate this on the basis of the ¹²C¹⁸O (3−2) transition (Fig. 4). Because of the low abundance of this carbon monoxide isotopologue, this line is less optically thick and thus traces deeper disk layers. However, its excitation temperatures of the lower energy levels are very similar to those of the ¹²C¹⁶O molecule. In contrast to the optically thick ¹²C¹⁶O line of the same model, the gap is clearly detectable.

The most massive disks (2.67 × 10^-1 M_⊙) in our sample do not show any signs of a gap in any of the considered transitions (see Table 3). For the reasons outlined above, for a specific molecule observation with the aim to reveal gaps, we recommend estimating the optical thickness of the disk (for example from its mass derived through continuum observations) in order to choose the appropriate molecule and (optically thin) transition.

Disk inclination

In contrast to the ideal face-on models, the velocity-channel maps of a disk with a small inclination of 10° is shown in Fig. 5. Negative velocities (top left corner) belong to parts of the disk moving toward the observer and areas with positive (bottom right) velocities are moving away from the observer. In this figure the main problem for detecting gaps in inclined disks, and thus the majority of potential targets, becomes obvious. Only fractions of the gap can be seen, but it does not appear as a symmetric ring as the gas covers a broad range of velocities (≫v_turb) along the gap throughout the individual velocity channels. This effect is illustrated in Fig. 6. The underlying disk model (Fig. 6, left) has an outer radius of 117 AU and the planet is located at 65 AU. Owing to the disk inclination of 10°, the line is significantly Doppler shifted. The velocity-channel map shows a cavity in the orbital region of the planet. However, this not an irrevocable sign of a gap. Owing to the projected velocity in the observer’s direction, the undisturbed model (see Fig. 6, right) also reveals a prominent local minimum in the same disk region. Thus, the velocity patterns effectively hide the gap and hinder unambiguous identification, especially for optically thick lines. A possible solution is to use velocity integrated intensity maps for gap detection instead (see Fig. 7). For low disk inclinations, these images do not deviate much from face-on inclined disks, i.e., for an inclination of 10° as considered in this study the flux deviation compared to the disk model seen face-on is less than 5%. The velocity integration has the advantage that the patterns disappear and the gap can be identified easily. However, the velocity dependency of the gap depth (see previous section) is lost. Consequently, the velocity-integrated maps are dominated by the line flux at the line center, which traces only the upper disk layers in the case of optically thick lines.

Central star

We now discuss the impact of using the Herbig Ae star, which is ~45 times more luminous, as the central star.

As the disk temperature is higher than the T Tauri case, the molecules are exited to higher energy states. The line radiation of the outer parts of the disk increases because a larger fraction of molecules is sufficiently exited to the relevant energy levels due to the warmer disk. In contrast to the situation in the outer disk regions, the higher temperature in the innermost region results in a depopulation of the lower energy states. Consequently, the low line transitions we consider in this study produce less flux in this region.

Second, the mass of the Herbig Ae star is about 3.5 times higher, resulting in a higher Keplerian velocity at a given radial distance from the central star. Thus, for the inclined disks the linewidth of the Doppler broadened line is much higher than in the T Tauri star case. Consequently, at a certain velocity channel (bandwidths of 50 m/s in this study) the parts of the disk with the appropriate velocities to contribute to this channel appear precise and well defined (see Fig. 8). Thus, assuming the same disks for both PMS stars, individual disk regions can be better studied in the case of the more massive central star, which increases the chances of gap detection. This effect is a direct result of the dependence of the linewidth on the kinetic temperature (cf. Eq. (B.8)).

Detectability of gaps in the ideal case

We describe our approach to identifying a gap and to deriving its depth in the velocity-channel maps.

We follow the procedure outlined by Ruge et al. (2013). We derive the contrast between the inner/outer edge of the gap and the gap itself by using radial cuts throughout the simulated maps (see Fig. 9). To reduce the probability of false detection (e.g., local minimum in the brightness distribution due to the inclination effect discussed in Sect. 4.1), 64 individual radial cuts are considered in the analysis. Only if a gap is indicated along more than 50% of successive radial cuts is it considered as detected. We successfully tested our algorithm with smooth undisturbed reference hydrodynamic models to verify that no substructure is misleadingly considered as a gap. However, because of the very general approach due to the huge parameter space, we emphasize that for specific observations better results may be achieved by adjusting the algorithm to the specific observation (e.g., a well-adapted compensation of background noise/sparse uv-coverage, which leads to asymmetric structures).

Overview charts

In this section, the feasibility of detecting gaps through observations of selected molecular lines for all considered disk configurations (see Table 4) is summarized. For this purpose, we derive the gap depth from the ideal velocity-channel maps for all five molecules and 32 assumed transitions using the algorithm described above. For every transition we derive the expected ALMA background noise of the reconstructed velocity-channel maps (σ) using the ALMA sensitivity calculator (see Sect. 4.2.1, Table 3). Owing to the atmospheric transmission, σ mostly rises with frequency resulting in a worse signal-to-noise ratio (S/N). To take care of this effect, we derive the gap depth and weight it by the expected sensitivity of ALMA at every line transition considered in this study. Thus, we obtain a criterion of how well individual transitions allow gap detection with ALMA. These results allow those transitions to be chosen that are best suited for the simulated ALMA observations in Sect. 4.2.

Figure 10 shows the results for an exemplary disk with an outer radius of 144 AU and a disk mass of 2.67 × 10^-3 M_⊙ around the Herbig Ae star. The gap depth measured in every velocity channel for the ¹²C¹⁶O (3−2), ¹²C¹⁸O (3−2), and CS (7−6) transition is shown. All transitions are optically thick at the line center. While the ¹²C¹⁶O (3−2) and CS (7−6) transitions fail to detect the gap at the line center, the optically thin ¹²C¹⁸O (3−2) transition allows for gap detection even at the central frequency. Nevertheless, the highest gap depth is achieved by the CS (7−6) transition at the line wings (at ±350 m/s).

To summarize our findings for all molecules, we search for the transition that provides the highest gap depth for given disk mass and disk size and present the result in one overview map. The resulting charts for both considered PMSs are shown in Fig. 11. The full overview charts showing the individual gap depth of all 32 transitions considered in this study can be found in Appendix A.

As already indicated in Sect. 4.1, the most promising molecule for identifying gaps in low-mass disks (~10^-4 − 10^-5 M_⊙) is ¹²C¹⁶O. The gap can be detected either at the line center for the lowest disk masses or at the line wings for the intermediate-mass disks. In most configurations the flux density and consequently the gap depth is higher even in the line wings compared to the gap depth which can be achieved with other molecules.

For the disks with intermediate masses (~10^-2 − 10^-3 M_⊙) of our sample gap detection becomes very challenging. In these cases ¹²C¹⁶O is optically thick for almost all considered disk configurations even at the line wings. Thus, transitions of less abundant molecules like HCO⁺, HCN, or CS are better suited for gap detection. However, because of the low abundance of these molecules the net flux is very low. Consequently, the gap depth that relates the flux density to the background noise is too low for gap detection in many cases. The most massive (~10^-1 M_⊙) disks in our sample do not show any reliable sign of a gap for any molecule transition considered. We find that the gaps of these disks can be traced best outside the line in the dust continuum, even with the low bandwidth considered in this study.

When comparing the molecular line emission from the disk around both PMSs, we find very similar results in terms of the probability of gap detection. The effects due to their different stellar mass and luminosity (see Sect. 4.1 for a discussion) lead to slightly different molecules, which offer the best probability of detecting a gap (e.g., because of different excitation states ¹²C¹⁶O allows for gap detection even for the disks with a mass of ~10^-4 M_⊙).

Therefore, our results demonstrate the feasibility of gap detection using molecular lines, but it is only efficient for low-mass protoplanetary disks or disks with huge inner cavities (e.g., transitional disks or circumbinary disks).

4.2. Simulated ALMA observations

4.2.1. ALMA simulation setup

We use the Common Astronomy Software Application package (CASA 4.2, Petry & CASA Development Team 2012) to simulate observations based on our ideal velocity-channel maps. We consider a total of six medium and extended configurations (baselines from 0.7 to 16 km) using 50 antennas, a total observing time of 3 h, a target distance of 140 pc, and the sky-coordinates of the “Butterfly Star” (α = 04^h33^m17^sδ = + 22°53′20′′, Wolf et al. 2003; Gräfe et al. 2013) which is representative of young stellar objects in the Taurus-Auriga star-forming region. The precipitable water vapor (PWV) in the atmosphere has a major influence. Therefore, we also include thermal noise in our simulations and use the recommended column densities for the line transitions we consider in our sample (see Table 3). First, we employ the simobserve task to generate realistic visibility measurements. The final simulated observations are obtained using the simanalyze task which uses the CLEAN-algorithm for image reconstruction (Högbom 1974). The expected noise levels σ are estimated with the ALMA sensitivity calculator³ because the measured noise in the simulated images strongly depend on the quality of the deconvolution⁴. The properties of the considered molecules and the expected sensitivity are summarized in Table 3. We consider the ALMA polarization configuration in dual-mode. With respect to the results from Sect. 4.1 we select those molecules that offer the highest probability of gap detection, i.e., ¹²C¹⁶O(3−2), ¹²C¹⁸O(3−2), CS (7−6), HCO⁺(3−2), and HCN (3−2). In addition, their corresponding transitions around 345 GHz (ALMA band 7, see Table 3), offer a reliable comparability of the different molecules and yield a good compromise between spatial resolution (higher frequencies) and better atmospheric transmission (lower frequencies)⁵. The ¹²C¹⁶O (4−3) transition is added to provide a second ¹²C¹⁶O transition offering higher spatial resolution at the cost of sensitivity.

4.2.2. Simulated line observations

In Figs. 12 and 13 a face-on and an inclined disk model are shown, respectively, as they could be observed with ALMA using configuration 14 (~1.6 km baseline) within three hours (corresponding ideal maps are shown in Figs. 2 and 5, respectively). More compact baselines result in higher flux densities and consequently higher S/Ns, while larger baselines allow gaps of compact disks to be traced at the cost of decreasing S/N. We find for the most extended disk configurations (R_in ~ 45 AU, R_out ~ 200 AU) in our model setup, maximum baselines of about 1 km are sufficient to identify the gap in these disks, but fail for the smaller disks.

For the compact disks in our parameter space, higher spatial resolution and thus longer baselines are required. In our sample the smallest disks have an inner/outer radius of 14 AU/63 AU, respectively. We find that unambiguous gap detection is feasible using intermediate baselines of about 1−3 km. Longer baselines are not required and only result in sparser uv-coverage and thus a lower S/N of the reconstructed maps. Furthermore, many of the smaller planetary signatures (inner and outer spiral arms and disk inhomogenities) are eliminated/smoothed out because of beam convolution and worse uv-coverage. However, for this particular model, and for most of the expanded models in our sample, the gap can be identified in both inclination cases using ALMA configuration 14.

For the same disk, the resulting HCO⁺ (4-3) velocity-channel map for ALMA configuration 14 is shown in Fig. 14. The HCO⁺ molecule is less abundant than ¹²C¹⁶O and consequently the line is less optically thick. Therefore, regions closer to the disk midplane can be observed and the gap appears more pronounced. In contrast, the line flux itself is about 35% weaker than the simulated ¹²C¹⁶O (3−2) line of the same model.

4.2.3. Gap detection results

Our main intention is to investigate the feasibility of tracing gaps in molecular lines with ALMA. To allow this for a large range of disk models, we define a general criterion for a quantitative measure of gap detection in the simulated ALMA maps. We apply the same algorithm as in the case of ideal maps (see Sect. 4.1). The resulting gap depth is given in units of σ to clarify the meaning of the gap detection. Only if the depth of the gap is greater than three times the noise level σ do we consider the gap as detected.

First, we demonstrate that the optical depth effect (see Sect. 4.1 for a discussion of the ideal case) has direct observational impact. As shown in the ideal case, the choice of the molecule for a given disk mass/size is most important for unambiguous gap detection. In Fig. 15 we present the gap depth for the same transitions as in Fig. 10, but now for a simulated ALMA observation utilizing configuration 14 (~1.6 km baseline). Both results are closely comparable as the gap depth in the ideal case shows the same trend as the simulated ALMA observations, i.e., ¹²C¹⁶O (3−2) cannot be used as a tracer for gaps, and ¹²C¹⁸O (3−2) allows for gap detection even at the line center, but the CS (7−6) transition traces the gap best and results in a gap depth of about 17σ at ±350 m/s.

Second, our results are summarized in the overview charts for the Herbig Ae (left) and T Tauri (right) star (Fig. 16). Compact ALMA configurations (9, 11) result in high S/Ns, but lack of spatial resolution. As a consequence, utilizing the most compact ALMA configuration nine gaps are only identifiable for the most extended models of our sample (e.g., disk models with R_out = 255 AU). If we consider longer baselines (higher spatial resolution) in the range of ~1.6−2.3 km we are able to identify gaps for more compact models with a σ> 3. Intermediate baselines (e.g., configuration 11−14, ~1.0−1.6 km) potentially allow gaps to be traced for the wide range of models considered in this study.

The comparison between the molecules considered shows that CS, HCN, and HCO⁺ allow the detection of gaps even for disks with intermediate masses (~10^-3 M_⊙), while ¹²C¹⁶O tends to fail. As already discussed in Sect. 3.2, owing to the optical thick ¹²C¹⁶O line only the surface of the disks can be observed, which is barely perturbed by the planet-disk interaction. In contrast, for lines of molecules with lower abundances the surface is still optically thin even for more massive disks. We find that especially CS and HCN are proper candidates to serve as a tracer for gap detection in the mass range from ~10^-3 M_⊙ to 10^-4 M_⊙.

For more massive disks (M_disk ≥ 10^-2 M_⊙) the ¹²C¹⁶O isotopologue ¹²C¹⁸O is an adequate candidate. However, owing to the resulting low line flux gaps these massive disks should be observed in the continuum because in the submm wavelength range provided by ALMA the disk is still optically thin in the gap region (see Ruge et al. 2013).

Finally, we find that because of the assumed observation time of 3 h, long-baseline configurations (20, 28) are not suited to identifying gaps because of the low intrinsic flux and inadequate uv-coverage.

5. Summary

Our intention was to explore how planet-induced structures, in particular planetary gaps, appear in high-angular resolution molecular line maps, and subsequently if ALMA would allow the detection of these structures. For this reason, we developed the line radiative transfer code called Mol3D. It allows the calculation of the level-populations with several N-LTE approximate methods and the production of ideal velocity-channel maps. First, we tested the new code against other established (molecular line) RT codes.

Based on 3D hydrodynamic simulations carried out with the PLUTO code, we simulated synthetic velocity-channel maps for five molecules (¹²C¹⁶O, ¹²C¹⁸O, HCO⁺, HCN, and CS) for a total of 32 rotational transitions. We summarize our findings:

• The patterns of the global (Keplerian) velocity field in thecorresponding channel maps of inclined disks maskplanet-induced structures. Gap identification can be ambiguousin these cases.

• Molecules that tend to become optically thick (like ¹²C¹⁶O) only trace the upper layers of the disk, which are rarely affected by the planet-disk interaction. In consequence, gap detection for massive disks (≳10^-2 M_⊙) is barely feasible with molecular lines.

• For lower mass disks, depending on the molecule and its optical properties (abundances, excitation state) the gap can be observed for a broad range of considered disk configurations, not necessarily at the central wavelength of this transition but at the line wings. Because of the higher optical depth of the line at its center, the τ = 1 surface traces only the upper unperturbed disk layers, while the τ = 1 surface of the line wings traces deeper and thus more strongly perturbed layers.

• With ALMA it is feasible to identify gaps using molecular lines for several disk configurations, if one respects the effects we outlined above (i.e., by choosing an appropriate optically thin molecule and corresponding transition). In particular, for low-mass disks (≲10^-3 M_⊙) the feasibility of identifying gaps using molecular lines is very promising. For this mass range molecular line observations fill the discovery space for gaps where continuum observations are supposed to fail (see Fig. 8 in the complementary study of Ruge et al. 2013).

• Our simulations are based on an assumed distance of 140 pc and a fixed sky coordinate (max. elevation ~45°) to give reliable results for the Taurus star-forming region. Thus, for many disks located closer or at more favorable positions on the sky the feasibility of detecting gaps is expected to be even higher.

We emphasize that the observation of molecules with different abundances allow for unambiguous identification of gaps and the physical and chemical composition of protoplanetary disks. Especially for low-mass disks where gap detection in the continuum is difficult owing to the low intrinsic flux, molecular line observations are a promising alternative with which to identify gaps.

Finally, we note that our study is based on HD simulations. In addition, Uribe et al. (2011) have remarked that magnetohydrodynamics (MHD) simulations result in more extended gaps owing to the presence of the weak magnetic field. Including intrinsic stellar radiation, the form of the gaps also changes because the gap walls can be directly radiated by the host star, resulting in deeper gaps. This effect is even higher for the outer wall (Jang-Condell & Turner 2012). Using smoothed particle hydrodynamic (SPH) simulations, Gonzalez et al. (2012) found that gaps are deeper for larger dust grains and bigger planets. In addition, we only considered a Jupiter-mass planet in this study. More massive planets would enlarge the gap width owing to the higher gravitational potential. If these effects are present it would therefore be even easier to observe gaps with ALMA than proposed within this study.

---

Acknowledgments

This presented work is funded through the DFG grant: WO 857/10-1. Special thanks go to Yaroslav Pavlyuchenkov for his comments and advice, which were very helpful in the early development phase of Mol3D. We also want to acknowledge Peter Scicluna and Tommy Breslein for language editing and Jan Philipp Ruge for the helpful advice, especially regarding the use of the CASA simulator.

References

Online material

Appendix A: Gap depth − complete overview charts

The full overview charts for all molecules considered in our study are summarized in Fig. A.1 for the T Tauri case and in Fig. A.1 for the Herbig Ae case. The derived σ-weighted gap depth is color coded for each considered disk configuration.

Appendix B: Mol3D

In this section we present a detailed derivation of the underlying methods, physics, and assumptions used by our new N-LTE line and dust continuum radiative transfer code Mol3D.

Appendix B.1: Radiative transfer equations

We start with the transfer equation in the following form, (B.1)with the monochromatic intensity I_ν, the emission factor j_ν(s), and the total absorption factor α_ν(s). The intensity from point p₁ to p₂ along its path s can be calculated using the integral form of Eq. (B.1), (B.2)where τ_ν is the optical depth between two points, p₁ and p₂, along the ray of sight defined as (B.3)

Appendix B.2: The gas-phase

For the gas, e.g., for a selected transition from lower level i to upper level j, the emission (j_ν) and absorption (α_ν) coefficients have the following form: We follow the ansatz of Rybicki & Hummer (1991) for the non-overlapping multi-level line treatment. For a molecule with N levels and a selected transition between lower i and upper j level, we consider the spontaneous downward transition rates A_ij, the Einstein coefficients B_ij for stimulated transition, and the collision rate C_ij. The collision rate can be calculated from the number density of the collision partner and the downward collision rate coefficient γ_ij, which is the Maxwellian average of the collision cross section σ measured in laboratory: (B.6)All molecule data we use in our software package (e.g., frequencies, transition rates A_ij, collision rate coefficients) are obtained from the Leiden Atomic and Molecular Database (LAMDA) (Schöier et al. 2005). As a consequence, Mol3D uses the same common input format for the molecule data as many other line radiative transfer programs available like RADEX (van der Tak et al. 2007) or URAN(IA) (Pavlyuchenkov & Shustov 2004). This approach offers the advantage of easily extending the code for several common molecules available in the LAMDA database (3 atomic and 33 molecular species at the time of writing). The dominant broadening effects are thermal and turbulent broadening. Hence, we assume a Gaussian line profile function, but in principle any other profile function can be included as well, (B.7)where ν_ij is the central frequency of the transition and c the speed of light. The parameter v_tot gives the total line broadening factor taking into account the assumed microturbulent velocity v_turb and the kinetic velocity v_kin due to the kinetic temperature T_kin: (B.8)In our code, the microturbulent velocity v_turb can either be set to a common value, ~0.1−0.2 km s^-1 for protoplanetary disks (Piétu et al. 2007; Hughes et al. 2011), or for example in the case of an outcome of (M)HD/MRI simulations, it can be provided as an input for every grid cell.

Finally, to describe the full line radiative transfer problem, a set of balance equations is needed: (B.9)In fact, Eq. (B.9) is not one single equation, but a whole set of linear equations, one for every energy level. In our code this linear matrix equation is solved with a simple Gaussian elimination algorithm. This set has to be solved locally (for every grid cell), but it has a global character due to its dependency on the line integrated mean intensity defined as (B.10)Equations (B.9) and (B.10) are directly coupled to Eq. (B.1).

To solve the line radiative transfer problem, one needs to start with estimated values for the level populations, solve the radiative transfer equation, calculate the mean intensity, and then solve Eq. (B.9) again to create a new set of level populations. Then one needs to iterate until the correct level populations have been found.

Appendix B.2.1: Level populations

The main problem in radiative line transfer is the calculation of the level populations. Equation (B.9) shows that the level populations are directly coupled to the mean intensity and therefore to the gas temperature, which is unknown for complex structures as in the case of protoplanetary disks.

We implemented dissimilar algorithms to calculate the level populations with three approximate methods, for example, local thermodynamical equilibrium LTE, FEP, and LVG. As Pavlyuchenkov et al. (2007) have shown for the case of protoplanetary disks, for the most common molecules LVG is a very accurate method of calculating the level populations and gives a very good ratio between computational efficiency and reliability of the results.

LTE level populations

In this assumption the level populations follow a Boltzmann distribution and the excitation temperature of the selected transition is assumed to be equal to the kinetic temperature T_exc = T_kin: (B.11)Here g_i and g_j are the statistical weights of the i − j line transition and ν_ij the corresponding central frequency.

LVG level populations

In the case of protoplanetary disks, the radial velocity gradients are usually much larger than the local thermal and microturbulent velocities (Weiß et al. 2005; Castro-Carrizo et al. 2007). This means that photons emitted at a certain disk region can only interact and consequently get absorbed locally. In this approach, the mean intensity J is approximated from the local source function S in combination with an external radiation field J_ext: The quantity β gives the probability of a photon beeing able to escaping the model. We note that in the optically thin case a photon can escape the model without interaction. For this case, β equals 1 and the LVG method would be equal to the full escape probability method. Therefore, the FEP method (also included in Mol3D) uses the same formalism as the LVG method, but with β always set to 1 (see also van der Tak et al. 2007). To calculate the local β, our program uses Eq. (B.14) introduced by Mihalas et al. (1978); de Jong et al. (1980) for accretion disks in general, (B.14)where τ is the effective optical depth of the observed line. We note that this formula is well suited in the case of proto-circumstellar disks. For other geometries and scenarios other approximations might give better results (see, e.g., de Jong et al. 1975 or Osterbrock & Ferland 2006). At this point an approximation of the optical line depth is made. As it is connected to the local velocity field V, microturbulence v_turb, and gas density, we use the following expression (see also Pavlyuchenkov et al. 2007): (B.15)

Appendix B.3: Benchmarks

We demonstrate the applicability of the code Mol3D. For this purpose we compare its results to those obtained with the line radiative transfer code URAN(IA) (Pavlyuchenkov & Shustov 2004) and the dust continuum radiative transfer code MC3D (Wolf et al. 1999; Wolf 2003).

Appendix B.3.1: URAN(IA)

URAN(IA) is a 2D radiative transfer code for molecular lines. It features a Monte Carlo algorithm for calculating the mean intensity and an accelerated lambda iteration (ALI) method for self-consistent molecular excitation calculations. This code has been extensively tested in 1D and 2D, for example against the RATRAN code (Hogerheijde & van der Tak 2000). It successfully passed all tests formulated in the benchmark paper for N-LTE molecular radiative transfer by van Zadelhoff et al. (2002). It has been used in several applications, for example modeling the starless core L1544 (Pavlyuchenkov & Shustov 2004), and in dedicated parameter space studies of molecular line formation of molecular line formation of prestellar cores (Pavlyuchenkov et al. 2008).

In this code, several approximate methods of calculating the level populations are included. In a comprehensive study, Pavlyuchenkov et al. (2007) have shown the reliability of these methods in protoplanetary disk scenarios.

In Sect. 2 we introduced the different methods of calculating level populations implemented in Mol3D. In Fig. B.1 we present the spectrum of the HCO⁺ (4−3) transition of a typical protoplanetary disk (M_disk = 0.07 M_⊙), comparing the different methods assuming a uniform N(HCO⁺)/N(H) ratio of 1 × 10^-8. For simplicity, we choose a flared disk model around a T Tauri star utilizing a temperature distribution with a gradient in radial and vertical direction and a power-law model for the density distribution adapted from Shakura & Sunyaev (1973): (B.16)The disk is orientated face-on. We assume pure Keplerian rotation and neglect dust re-emission. All spectra are normalized to the maximum value of the LTE solution. To allow comparison, the resulting spectra obtained with the URAN(IA) code are also included.

URAN(IA) and Mol3D produce comparable HCO⁺ (4−3) spectra with differences of less than about 0.5%. The most realistic solution (black line) is obtained using the accelerated Monte Carlo (ART), to classify the approximation methods. As discussed in greater detail in Appendix B.2.1, the LTE method overestimates the net flux and the FEP method underestimates it. The fluxes in the line obtained with the LVG method are in between those obtained with the other methods and are closely comparable to the ART method. This result has also been found by Pavlyuchenkov et al. (2007) and we refer the interested reader to their study.

Appendix B.3.2: MC3D

MC3D is a 3D continuum radiative transfer code that uses a Monte Carlo method to calculate self-consistent dust temperature distributions, continuum re-emission/scattering maps, and SEDs. Among a wide variety of radiative transfer studies it has been applied extensively for the analysis of multiwavelength high-angular-resolution observations of circumstellar disks (e.g., Schegerer et al. 2008; Sauter & Wolf 2011; Madlener et al. 2012; Gräfe & Wolf 2013). MC3D was successfully tested against other continuum RT codes (e.g., Pascucci et al. 2004).

In this section, we compare Mol3D with MC3D for a typical protoplanetary disk scenario. The density distribution is described by Eq. (B.16), assuming a total disk mass of 0.04  M_⊙. The temperature structure is calculated self-consistently with

each code using their Monte Carlo scheme. Both codes use a similar approach based on the assumption of local thermal equilibrium and immediate temperature correction as proposed by Bjorkman & Wood (2001).

The resulting midplane temperature around a PMS T Tauri star is shown in Fig. B.2. The calculation of the temperature in the dense midplane is one of the most crucial RT problems, because these disk regions can hardly be reached by stellar photons directly. Thus, it is mostly heated indirectly by thermal re-emission radiation from the innermost disk regions and upper disk layers. Both codes produce very smooth and similar temperature distributions. The statistical nature of the Monte Carlo method results in deviations of about 10% for the innermost disk parts and less than 2% for the outer regions, which is mainly due to the probability that the photon packages can reach these regions.

We also compared the dust re-emission/scattered light maps and SEDs and find that both codes produce comparable results with maximum deviations of about 10%.

All Tables

All Figures