© ESO, 2013

1. Introduction

Disks rich in dust and gas around recently-formed stars are important because they harbor either recently-formed planets or possibly those still in the process of forming. The structure of these young protoplanetary disks has been the subject of intense study over a wide range of wavelengths (see Williams & Cieza 2011 for a recent review). Their sizes range up to a few hundred AU with Keplerian rotation velocities and temperatures up to a few tens of K in the disk midplane several AUs from the central star. This means that millimeter/submillimeter telescopes are well suited to study their molecular and dust components. Single-dish telescopes generally cannot resolve the disks, but have revealed the characteristic double-peaked line profiles from CO and several other species. Interferometers at millimeter wavelengths have provided images a few resolution elements across, but have not yet provided highly detailed images.

At a distance of 122 pc, HD 163296 is an isolated A2Ve star with an estimated age of about 5 Myr (Montesinos et al. 2009). It is one of the most thoroughly studied protoplanetary disks, and was one of the first to be resolved with millimeter interferometry (Mannings & Sargent 1997). Submillimeter interferometry observations with 2″ resolution indicates an outer radius of 550 AU in CO, with an inclination of 45° (Isella et al. 2007). The brightness of this source in millimeter molecular lines has made HD 163296 an excellent laboratory for the comparison of disk models, inspiring studies of radial and vertical temperature and molecular abundances (Qi et al. 2011; Tilling et al. 2012; Akiyama et al. 2012). Determining unique solutions to the underlying disk structure clearly benefits from having high-angular resolution, good image fidelity, and high-sensitivity. Detailed studies of dust-and-gas structures in protoplanetary disks are now possible with the spatial resolution and sensitivity provided by ALMA.

In this work we present a detailed study of the disk structure surrounding HD 163296 using the best images to date provided by ALMA band 7 data in continuum and spectral line.

2. Observations description

Observations were performed on 2012 June 9, 11, 22, and July 6 at Band 7 as part of the ALMA science verification program 2011.0.000010.SV. The array was in a configuration with projected baselines length between ~16 to ~400 m, sensitive to maximum angular scales of ~7′′ and providing a synthesized beam of 0.52″ × 0.38″ at PA ~ 82°. The field of view was ~18′′. A total of five data sets were collected, using between 16 and 19 antennas of 12 m diameter and accounting for 3.9 h of total integration time including overheads and calibration (2.3 h on HD 163296). Weather conditions were good and stable, with an average precipitable water vapor of 0.8 mm. The system temperature varied from 100 to 300 K.

The correlator was set up to four spectral windows in dual polarization mode, centered at 345.796 GHz (CO(3–2)), 346.998 GHz (H¹³CO(4–3)), 356.734 GHz (HCO⁺(4–3)), and 360.170 GHz (DCO⁺(5–4)). The effective bandwidths used were 468.75, 937.50, 468.75, and 117.19 MHz, respectively, at velocity resolutions of ~0.21, 0.42, 0.21, and 0.05 km s^-1 after Hanning smoothing. In this work we present results from the observations in the CO(3–2) emission line and the continuum at ~850  μm.

The ALMA calibration includes simultaneous observations of the 183   GHz water line with water vapor radiometers that measure the water column in the antenna beam, later used to reduce the atmospheric phase noise. Amplitude calibration was done using Neptune, and quasars J1924− 292 and J1733− 130 were used to calibrate the bandpass and the complex gain fluctuations, respectively. Data reduction was performed using the version 3.4 of the Common Astronomy Software Applications package (CASA). We applied self-calibration using the continuum and we used the task CLEAN to image the self-calibrated visibilities. The continuum image was produced by combining all of the line-free channels. Briggs weighting was used in both continuum and line images. The achieved rms was 0.5 mJy beam^-1 for the continuum and 14 mJy beam^-1 for each CO channel map.

3. Results and modeling

3.1. Continuum emission at 850 μm

Continuum emission at 850 μm is detected centered at the position R.A.(J2000) = 17^h56^m21285, Dec(J2000) = −21°57′22″368. The dusty disk is well resolved with no suggestion of gaps or holes at radius > 25 AU. The major axis has a projected diameter of 3.9′′ (measured at the 3σ level), which corresponds to an outer dust disk radius of R_out ≃ 240 AU at the adopted distance of 122  pc (van den Ancker et al. 1998). This value is ~20% larger than the one reported by Isella et al. (2007) and can be explained by the higher sensitivity and resolution of the ALMA data, unveiling the fainter outer edge of the dusty disk (see Fig. 1). The inclination angle is i = 45°, derived from the shape of the isophote contours, and it is in agreement with Isella et al. (2007). The major axis position angle is 137°. The flux density integrated over the whole disk structure above 3σ is S_850  μm = 1.82 ± 0.09 Jy, similar to the values reported by Isella et al. (2007) and Qi et al. (2011) using Submillimeter Array (SMA) observations at similar wavelengths.

The continuum surface brightness profile was calculated by azimuthally averaging the emission from concentric elliptical annuli from the central star. It can be fitted by a two-component power law, with a characteristic radius (R_C, outside of which the brightness profile drops toward zero) equal to 125 ± 5 AU from the central star, changing from a power-law slope of 1.03 ± 0.04 to a very steep decline of 4.7 ± 0.2 (see Fig. 2). A minimized χ² fitting was used to determine the best value of the power-law slopes.

Fig. 1 Colors represent the CO(3–2) integrated emission over all velocity channels (first-order moment, 5σ cut). Contours show the continuum emission at 850 μm with levels representing 5, 10, 20, 40, 80, 160, and 320 times the rms of the continuum map (0.5 mJy  beam-1).

3.2. CO(3−2) spectral line emission

The CO data show a disk in Keplerian rotation (see Fig. 1). Considering contours above 3σ, the outer radius extends to R_out ≃ 575 AU, in agreement with Isella et al. (2007). The integrated intensity averaged over the whole emission area is 100.30 ± 0.13 Jy km s^-1, in agreement with the measurement reported by Qi et al. (2011). The inclination and the PA of the gaseous disk (measured at the 3σ contour level) are 38° and 138°, respectively. We note that the inclination value is different from the one using the continuum because of optical thickness and flaring combined with high-spatial resolution effects (see next section). For our analysis we adopt the values calculated from the continuum that are less affected by the flared structure of the disk.

The same fitting as done for the continuum in Sect. 3.1 was applied to the line data. In this case the surface brightness profile of the CO(3–2) can also be approximated by a power law with a break in the surface brightness at 125 ± 5 AU (similar to that of the continuum). At radii 125–240 AU, the CO surface brightness drops with a power law of 1.0 ± 0.1 and beyond 240 AU, the CO drops rapidly with a power law of 2.5 ± 0.1 (considerably less steeply than the continuum; see Fig. 2).

Fig. 2 Surface brightness profiles for the continuum (red crosses) and for the CO(3–2) (light blue crosses) with 1σ error bars. Dashed and dash-dotted lines represent two different model fits required for the continuum (with γ = 0.1) and for the spectral line (γ = 0.9) profiles respectively. Green dotted line marks Rc = 125 AU.

3.3. Modeling

To extend beyond the simple power-law fitting, the three-dimensional Monte Carlo radiative transfer code MCFOST was used to model the data presented in this work (see Pinte et al. 2006, 2009 for a more detailed description). The initial disk model we use is based on Tilling et al. (2012), in particular their favored model number 3 (see Table 6 in that paper). The same photospheric parameters were used, namely T_eff = 9250 K, M_∗ = 2.47  M_⊙, L_∗ = 37.7  L_⊙, and a slight UV excess, L_UV/L_∗ = 0.097. The model contains an exponentially tapered-edge and provided an adequate fit to the broadband spectral energy distribution (SED), the Herschel PACS lines observed by the Key program GASPS (Gas in Protoplanetary Systems; Dent et al. 2013), as well as submillimeter interferometry (continuum and CO data, from Isella et al. 2007).

The temperature structure in the disk is calculated by considering the dust opacity only, assuming astronomical silicates (Draine & Lee 1984) with a power-law size distribution having a slope of − 3.5 and a maximum radius of 1 mm. Each cell in the computation domain has its own density, grain properties, and opacity constructed following the global disk parameters. The calculated dust temperature profile (T_dust) decreases outward and there is a vertical temperature gradient for the dust that also depends on radius, with the midplane being cooler than the disk surface.

To calculate the CO(3–2) channel maps and surface brightness distribution MCFOST assumes a constant gas-to-dust mass ratio of 100 throughout the disk (both radially and vertically). We adopted a standard CO abundance with respect to H₂ (10^-4), set constant through the disk where T_dust > 20 K and equal to zero where T_dust < 20 K to mimic the effect of CO freeze out (see Sect. 4.3). The level populations are calculated assuming LTE and T_gas(r,z) = T_dust(r,z) for each grid cell. The radial and vertical temperature profiles and the radiation field estimated by the Monte Carlo simulation are used to calculate level populations for the CO molecule and to produce the SED, continuum images, and line emission surface brightness profiles, as well as kinematics with a ray-tracing method. The kinematics are calculated assuming the disk is in pure Keplerian rotation.

4. Discussion

4.1. Dust and gas surface brightness distributions

In order to fit the brightness radial profiles observed in the continuum at 850 μm and in the CO(3–2) emission line, we consider a tapered-edge model for the surface density distribution (see Andrews et al. 2009), (1)where R_c is the characteristic radius and γ is the index of the surface density gradient. In our data, CO is detected over a radius more than twice the dust continuum radius. This feature, for which the tapered-edge disk model has provided a solution for previous studies at lower spatial resolution and sensitivity, was produced because the CO line opacity is very much larger than the dust continuum opacity at 850 μm, such that the CO gas remains optically thick and detectable over a much larger radius.

4.1.1. The outer disk, outside of R_c

The initial disk model was based on the best model found by Tilling et al. (2012), which adequately fits to the broadband SED, the [OI] 63 μm Herschel PACS line, as well as previous submillimeter interferometric observations in CO(3–2), CO(2–1), and ¹³CO(1–0). To perform the fitting to the radial profiles of the CO(3–2) and the continuum, the R_c and the power-law slope derived from the simple initial fitting described in Sect. 3 were used as a starting point. The previous value of γ = 1, reported in the literature by Hughes et al. (2008) as the value that simultaneously fits the continuum and the CO(3–2) emission as observed with the SMA, was also considered as a reference for the fitting. A range of γ values from 0 to 2 were then explored to reproduce the steep decline in continuum and the shallower CO(3–2) surface brightness profile (see Table 1 for a description on the parameter space used).

Good model values for R_c and γ were found separately for the continuum and for the line surface brightness data. No minimization scheme was used, but we note that the continuum model and the line model (outside of R_c) match all the data points within the error bars (see discussion below for the inner 100 AU in the CO line). A posteriori, we verify that both models use the same value for R_c however, the line and continuum require very different values for γ; after the best γ value was found, the total flux density and the size of the emission provided by the model were compared with the values obtained observationally for both continuum and line. A dust mass of 7 × 10^-4  M_⊙ (Tilling et al. 2012) was assumed and a comparison of the synthetic SED with the observations was also performed for consistency.

As a result, the best match for the very steep decline of the continuum surface brightness outside R_c requires very small values of γ (<0.1; see Fig. 2), reproducing accurately the total integrated continuum flux (5% larger than the observed integrated flux), the R_c, and the radial profile of the surface brightness. But such a small gamma cannot fit the integrated CO(3–2) radial profile and the model line emission would only be detectable out to r ~ 400 AU, while it is observed out to 575 AU. Using γ = 0.8 − 0.9 provides a better match to the CO(3–2) surface brightness profile data, as well as for integrated CO flux (1% higher than the observed value), and R_c (see Fig. 2). However, a value of γ as high as this would result in a continuum disk whose surface brightness decline is too shallow. Values of γ near 1.0 have been reported in previous studies at lower resolution to fit well both the spectral line emission and the continuum (e.g., Hughes et al. 2008). Those values are ruled out for the continuum by the higher angular resolution data reported here.

Fig. 3 Channel maps of the CO(3–2) emission in HD 163296 from 3.19 to 8.23 km  s-1. The spectral resolution is 0.21 km  s-1 and the rms per channel is 14 mJy  beam-1.

4.1.2. The central region of the disk, inside of R_c

For the continuum data, the same model with small γ matches the observations nicely and no modifications are required. However, there is an extra component for the CO line emission at the center of the disk that is visible in Fig. 2. It also shows as a break in the CO(3–2) surface brightness profile just inside of R_c. To investigate the origin of this extra emission, we considered further the results of Tilling et al. (2012), which was the inspiration for our model. They showed (see their Fig. 7, right panel) that there is a layer of hot gas at the surface of the disk, significantly hotter than the dust, for radii smaller than 50–80 AU. The MCFOST code does not calculate the gas temperature, but we can easily verify that the model CO surface brightness profile is a good match to the data when using T_gas = 1.5 T_dust for a radius <80 AU, and T_gas = T_dust elsewhere. This is the model shown in Fig. 2, where the resulting change in the brightness profile, after convolution by the synthesized beam, is shown to match the observed profile. In this slightly modified model, the same value of gamma =0.9 is used for the CO gas disk, and the extra emission comes only from the increase in gas temperature at small radii, with no impact on the surface brightness profile of the outer disk.

Fig. 4 CO(3–2) representative channel emission maps in HD 163296, primary beam corrected. The spectral resolution is 0.21 km s-1 and the rms per channel is 14 mJy beam-1. The synthesized beam size is represented in the lower-right in each panel. The white solid and dotted lines represent the front and back surfaces of the gas disk, respectively. The insert at the bottom-left of each panel shows a zoomed-down model at the corresponding velocity (i.e., the size is the same, the scale is different).

4.2. Origin of the different surface brightness profiles

Coronagraphic imaging with the Hubble Space Telescope toward HD 163296 revealed a scattered light disk with an outer radius that extends to at least 450 AU (Grady et al. 2000). This suggests submicron dust particles responsible for scattered light remain coupled to the gas, while larger particles (~100 μm or larger) responsible for the 850 μm emission are concentrated in a smaller radius closer to the central star. One plausible explanation for the clear difference in the radial distributions of millimeter-sized dust grains and CO gas is a combination of grain growth and inward migration. It is known that grain growth occurs in this disk, as reported by Isella et al. (2007) based on the slope of the dust opacity law in the interval 0.87 − 7 mm. Related to the inward migration, models predict the gas drag to be more efficient in intermediate sized dust particles, with small grains remaining coupled to the gas (as seen in scattered light) and boulders following marginally perturbed Keplerian orbits (e.g., Barrière-Fouchet et al. 2005). Models predict a sharp cutoff of the submillimeter continuum emission when both grain growth and inward migration are considered (e.g., Laibe et al. 2008), in agreement with our observations. A similar case of a dusty compact disk with a sharp outer edge smaller than the CO gaseous disk has been reported by Andrews et al. (2012) for TW Hya.

Radial dust migration, expected to occur in disks massively, is an important ingredient for the evolution of disks and the formation of planets. Nevertheless it remains a phenomenon that is poorly constrained by observations. More efforts are needed to identify new sources like TW Hya and HD 163296 to better understand and characterize the physical processes leading to these more compact continuum disks at submillimeter wavelengths.

4.3. The vertical structure of HD 163296

The CO(3–2) channel maps unveil clearly for the first time at submillimeter frequencies the vertical structure of a flared gaseous disk (see Fig. 3). In Fig. 4 we have plotted four representative channel maps that show bright CO emission from the front disk surface, as well as similar but fainter and apparently rotated emission from the rear disk surface. This two-layer effect was predicted theoretically by Semenov et al. (2008) and it arises because the disk is tilted, the CO(3–2) is optically thick, and its emission is close to the disk surface where the flaring effect is more prominent.

To reproduce the spatial and kinematic pattern observed in the CO channel maps, the CO gas model derived in the previous section was used with γ = 0.9 and without the extra heated CO layer in the center. Because there is an intrinsic vertical temperature gradient calculated in the disk (with the midplane being colder), two-layer channel maps are produced naturally with lower intensity CO emission coming from the midplane compared to the surfaces. However, without further changes to the model, these two-layer channel maps do not produce enough contrast, i.e., enough difference in surface brightness between the midplane and the surface layers, even when convolved by the proper beam. Furthermore, the two surface layers appeared slightly too close to one another. Both effects are suggesting the need for a geometrically thicker disk and/or a larger vertical temperature gradient.

In order to improve the match between modeled and observed channel maps, we fine-tuned the CO model slightly by exploring a set of values of CO freeze-out temperature, non-thermal turbulence velocity, scale height, and flaring power. A summary of the parameter space explored is shown in Table 1.

First, we analyzed the effect of the variation of the freeze-out temperature in the patterns in the channel maps. The contrast between the warmer upper layers and the colder midplane can be increased by directly removing CO gas from the midplane. This can be done, for example, when the temperature is low enough that CO ice is produced, removing CO from the gas phase. To mimic the effect of freeze out in the midplane, we set the CO gas abundance to zero wherever T_gas = T_dust is below a critical value. Based on the initial value of 19 K adopted by Qi et al. (2011), a range of higher and lower values were explored. In particular a range of freeze-out temperatures from 10 to 40 K were examined. The best match is obtained for 20 K, similar to the value obtained by Qi et al. (2011). For 20 K, the layer where CO is removed in our model has a vertical extent of 15 AU (at 200 AU radius), and it affects mostly the outer disk midplane. This provides direct evidence for CO freeze out close to the disk midplane in HD 163296 (see also Mathews et al. 2013 in press). Freezing out CO gas improves the match between the model and the observations.

The effect of turbulence was then estimated by varying the non-thermal velocity component of the line width. Values similar to those estimated in other disks (e.g., Hughes et al. 2011) were considered. In particular, values between 0 to 0.3 km  s^-1 were explored and compared with the spatial extent of the emission in a given channel. The best match was obtained for 0.1 km  s^-1. This is slightly less than the value of 0.3 km  s^-1 suggested for the upper layers of the disk in HD 163296 by Hughes et al. (2011) for the CO(3−2) line, although a significant uncertainty is attached to this number. In our calculations, values higher than 0.2 km  s^-1 smear out the emission patterns too much in each channel.

The apparent separation between the two layers of CO emission in the channel maps is a function of the system’s inclination and geometry (thickness) of the disk. The inclination is well known, as derived from the observations. The geometry of the disk in our model is defined in a large part by the reference scale height (H₀) and the flaring exponent. Starting from the initial scale height and flaring index reported by Tilling et al. (2012), a range of higher and lower values of the scale height (between 0.03 and 0.1 at 1 AU) and flaring index (from 1.00 to 1.21) were considered (see Table 1). A good match is obtained for a scale height of 0.07 AU (at 1 AU radius) and a flaring power of 1.12. These values correspond to a disk that is slightly geometrically thicker than the disk of Tilling et al. (2012).

This set of parameters provides a combination of values that fit the observed CO(3−2) channel maps reasonably well and that are close to the values reported in previous studies. The final model also matches the CO radial surface brightness profiles adequately. Nevertheless, when the parameters of the refined model are used to fit the SED (i.e., a thicker disk in both the dust and gas), it is degraded in the mid- and far-infrared. This can be seen in Fig. 5 where 30% error bars are represented in the photometry measurements from mid- to far-infrared that cannot account for the discrepancy observed between model and observations. Thus, while the CO(3–2) channel maps require a thicker gas disk, the mid- and far-infrared part of the SED would require a flatter dust disk. A possible explanation of this discrepancy is that a significant amount of dust settling has occurred in the disk or, alternatively, that the disk is somehow stratified and there is a gas layer devoid of dust present on top of the dust disk a few scale heights above the mid-plane (the CO(3–2) τ = 1 surface is located at Z = 33 AU at R = 100 AU). In order to investigate further these phenomena, theoretical models that account for the effects of dust settling and radial migration would be needed.

Fig. 5 Spectral energy distribution of HD163296. Photometry values are taken from Tilling et al. (2012) and references therein. WISE data in all four bands are also included (Cutri et al. 2012). The blue line represents the fit based on the best model from Tilling et al. (2012). The dashed blue line represents our refined model that matches the CO(3−2) channel maps more accurately, with increased scale height and flaring. For better comparison ±30% error bars are plotted in red for the photometric measurements between 4.6 and 160 μm.

5. Conclusions

We have presented ALMA observations in CO(3–2) and continuum at 850 μm. We used these data in combination with MCFOST models to refine the geometrical, physical, and chemical properties of the disk surrounding HD 163296. The main conclusions are:

• For the first time at submillimeter frequencies, a detailed disk gas structure is shown, where the front and back surfaces of the gas disk are resolved.

• The high spatial resolution and sensitivity of the ALMA data presented in this work required refinements to the model previously presented by Tilling et al. (2012) that inspired our model. As a result, the CO(3–2) channel maps require a thicker gas disk in order to fit the observations, while a flatter dust disk is needed to reproduce the SED in the mid- to far-infrared wavelengths. We conclude that this could be evidence that the CO(3–2) emission comes from gas located farther above the disk midplane than the dust or, alternatively, an indication of dust settling in the disk. We confirm that the CO freeze out (at ~20 K) improves the fit to the CO data as well.

• We found a clear and strong difference in the radial surface brightness distribution of CO and submillimeter dust. A tapered-edge model is unable to fit both sets of data simultaneously with the same density profile for both. The critical radius of the models are similar, but the dust disk requires a much steeper cutoff and sharp outer edge. We conclude that the adopted tapered edge prescription, valid for a number of disks observed at lower resolution, must be modified to account for the effects seen with the higher sensitivity and angular resolution now available with ALMA. We propose a combination of grain growth and inward migration as a plausible explanation of this discrepancy.

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.000010.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. We are grateful to Dr. Antonella Natta who provided very useful comments. I.d.G. acknowledges the Spanish MINECO grant AYA2011-30228-C03-01 (co-funded with FEDER fund). F.M.e. and A.H. acknowledge financial support provided by the Milenium Nucleus P10-022-F, funded by the Chilean Government, and by the EU FP7-2011 programme, under Grant Agreement 284405. C.P. acknowledges funding from the European Commission’s 7th Framework Program (contract PERG06-GA-2009-256513) and from Agence Nationale pour la Recherche (ANR) of France under contract ANR-2010-JCJC-0504-01.

References

All Tables

All Figures

	Fig. 1 Colors represent the CO(3–2) integrated emission over all velocity channels (first-order moment, 5σ cut). Contours show the continuum emission at 850 μm with levels representing 5, 10, 20, 40, 80, 160, and 320 times the rms of the continuum map (0.5 mJy  beam-1).
In the text

	Fig. 2 Surface brightness profiles for the continuum (red crosses) and for the CO(3–2) (light blue crosses) with 1σ error bars. Dashed and dash-dotted lines represent two different model fits required for the continuum (with γ = 0.1) and for the spectral line (γ = 0.9) profiles respectively. Green dotted line marks Rc = 125 AU.
In the text

	Fig. 3 Channel maps of the CO(3–2) emission in HD 163296 from 3.19 to 8.23 km  s-1. The spectral resolution is 0.21 km  s-1 and the rms per channel is 14 mJy  beam-1.
In the text

Fig. 4 CO(3–2) representative channel emission maps in HD 163296, primary beam corrected. The spectral resolution is 0.21 km s-1 and the rms per channel is 14 mJy beam-1. The synthesized beam size is represented in the lower-right in each panel. The white solid and dotted lines represent the front and back surfaces of the gas disk, respectively. The insert at the bottom-left of each panel shows a zoomed-down model at the corresponding velocity (i.e., the size is the same, the scale is different).

In the text

Fig. 5 Spectral energy distribution of HD163296. Photometry values are taken from Tilling et al. (2012) and references therein. WISE data in all four bands are also included (Cutri et al. 2012). The blue line represents the fit based on the best model from Tilling et al. (2012). The dashed blue line represents our refined model that matches the CO(3−2) channel maps more accurately, with increased scale height and flaring. For better comparison ±30% error bars are plotted in red for the photometric measurements between 4.6 and 160 μm.

In the text