A SPATIALLY RESOLVED VERTICAL TEMPERATURE GRADIENT IN THE HD 163296 DISK

Spectral line emission from molecules can be a powerful tool for studying the structure of protoplanetary disks. Abundant molecules in disks, such as CO or CN, trace the spatial distribution of gas temperatures and the radial extent of the molecular hydrogen that composes most of the disk mass (Koerner et al. 1993; Guilloteau et al. 2013). Analysis of multiple lines can determine the vertical temperature structure of the disk (Dartois et al. 2003) or indicate freezeout of the gas phase molecules onto dust grains in the cold midplane (Qi et al. 2011). Observations of various species are useful for studying the abundance distributions of molecules and processes like ionization (Öberg et al. 2011), grain surface reactions (Dutrey et al. 2011), and fractionation (Öberg et al. 2012). Furthermore, these spectral lines probe the bulk motions, or kinematics, of the gas (Beckwith & Sargent 1993). Observations from interferometers are particularly useful since the emission is resolved both spatially and spectrally and can be used to derive a dynamical mass of the young star (Guilloteau & Dutrey 1998; Simon et al. 2000), probe the structure of the disk (Dutrey et al. 2008; Panić et al. 2010), or detect non-thermal line broadening from turbulence (Hughes et al. 2011; Guilloteau et al. 2012).

The Herbig Ae star HD 163296 (spectral type A1) hosts a protoplanetary disk with bright continuum emission (Allen & Swings 1976; Mannings & Sargent 1997; Natta et al. 2004) and a rich molecular spectrum (Qi 2001; Thi et al. 2001, 2004; Qi et al. 2013). It is an isolated system, not known to be associated with any star forming region or young moving cluster (Finkenzeller & Mundt 1984). Observations of CO lines show that the disk exhibits Keplerian rotation (Mannings & Sargent 1997; Isella et al. 2007; Hughes et al. 2008) and is remarkably large, extending past a radius of ∼500 AU (Isella et al. 2007). The disk is seen in scattered light (Grady et al. 2000) and is associated with an asymmetric Herbig–Haro outflow (Devine et al. 2000) and molecular disk wind (Klaassen et al. 2013). The system is young (∼5 Myr; Natta et al. 2004) and its Hipparcos parallax indicates it is nearby ($d= 122^{+17}_{-13}$ pc; van den Ancker et al. 1998). Qi et al. (2011) developed a density and temperature structure for this disk that was consistent with Submillimeter Array (SMA) observations of both the dust continuum and multiple CO emission lines. Including an abundance distribution that accounts for both freezeout and photodissociation, their analysis suggested that the disk has a midplane CO snow line at r ∼ 155 AU. An ancillary result of their model is that the disk is colder in the midplane than in the atmosphere where the higher J emission lines are produced. Evidence for a vertical temperature gradient in this disk was independently presented by Akiyama et al. (2011).

We present an analysis of the spectral line emission from the main CO isotopologues observed by the Atacama Large Millimeter/sub-millimeter Array (ALMA) toward the HD 163296 disk in Bands 6 and 7. We first use a toy model comprised of a rotating double cone to guide our intuition of this data. Leveraging the high spatial and spectral resolution of these observations, we then demonstrate that the vertical temperature gradient in the disk is directly resolved in the ¹²CO J = 3–2 maps alone. We develop a disk structure model that is roughly self-consistent with those data and the J = 2–1 transitions of ¹²CO, ¹³CO, and C¹⁸O. In addition, we explore the velocity structure of the disk and show that the excellent spatial and spectral resolution of the sensitive Band 7 observations can distinguish between various classes of disk models which have bulk gas velocities that deviate at the ≳ 5% level. We describe the observations in Section 2 and show the results in Section 3. Our analysis follows in Section 4, where we present the vertical temperature gradient (Section 4.2) and explore sub-Keplerian motions due to the vertical geometry of the disk and its radial pressure gradient (Section 4.3). We conclude with a brief discussion in Section 5 and a summary in Section 6.

HD 163296 ($\alpha = 17^{\rm h} 56^{\rm m} 21\buildrel{\mathrm{s}}\over{.}287$, δ = −21°57'2239, J2000) was observed by ALMA in Band 6 and Band 7 as part of its commissioning and science verification (SV) program. The raw and calibrated visibility data were publicly released through the science portal, and were accompanied by a set of detailed calibration and imaging scripts provided by the ALMA SV team.³ We started with the calibrated measurement set and used the CASA software package (v3.4; McMullin et al. 2007) to produce self-calibrated and continuum subtracted spectral visibilities for the ¹²CO, ¹³CO, and C¹⁸O J = 2–1 lines at 230.538, 220.399, and 219.560 GHz in Band 6 and ¹²CO J = 3–2 at 345.796 GHz in Band 7.

The Band 6 observations were taken using 24 of the ALMA 12 m antennas on 2012 June 9, June 23, and July 7 for a total on-source time of 84 minutes (including latency). The baselines spanned a range of 20–400 m. The correlator was configured to simultaneously observe four spectral windows (SpWs), with two in each sideband: SpWID 1 included the ¹³CO and C¹⁸O J = 2–1 transitions, while SpWID 2 covered the ¹²CO J = 2–1 line. SpWIDs 0 and 3 observed line-free continuum (but see Qi et al. 2013) and were centered at 217 and 234 GHz. The flux calibrator for each of the nights was Juno, Neptune, and Mars, respectively, and the bandpass solution was generated from observations of the quasar J1924-292. Integrations on the science target were interleaved with ∼2 minute long observations of the nearby quasar J1733-130 for phase calibration. For the July 7 observations, the flux density of J1733-130 was determined to be 1.55 Jy, using a bootstrap from the flux calibrator, Mars. Similarly, HD 163296 was observed in Band 7 on 2012 June 9, June 11, June 22, and July 6 with the same antenna configuration. Juno (June 9 only) and Neptune were observed as primary flux calibrators. When Neptune was observed, the flux density of J1733-130 was bootstrapped for each individual SpW. An additional bootstrap from these observations was used to determine the flux density of J1733-130 in each SpW observed on June 9. The total on-source time was 140 minutes. These observations used the same bandpass and phase calibration strategy as the Band 6 observations. The correlator was set up to observe four SpWs, with two in each sideband: SpWIDs 0 and 3 were centered at 360 and 347 GHz for continuum, HCO⁺ J = 4–3 at 356.734 GHz was detected in SpWID 1, and the ¹²CO J = 3–2 line was covered by SpWID 2.⁴

Figure 1 shows channel maps of the ¹²CO J = 3–2 line along with the J = 2–1 lines from ¹²CO, ¹³CO, and C¹⁸O. Before imaging, the visibilities have been regridded onto the same 0.34 km s⁻¹ wide channels; only the central 15 channels are shown, where the emission pattern is clearly resolved for all lines. The central 49 channels of the ¹²CO J = 3–2 line imaged at the native Δv = 0.11 km s⁻¹ channel spacing is also shown in Figure 2. We generated these synthesized images using the CASA implementation of the CLEAN algorithm (Högbom 1974) with natural weighting. Table 1 summarizes the noise characteristics, synthesized beams, native channel widths, and integrated intensities for each line. The reported line intensities agree with previous measurements by the SMA (Qi et al. 2011) within the 10% systematic uncertainty estimated for the absolute calibration of the visibility amplitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Emission Line Results

Parameters	12CO(3–2)	12CO(2–1)	13CO(2–1)	C18O(2–1)
Channel widtha (km s−1)	0.11	0.32	0.33	0.33
Beam	065 × 042	081 × 066	087 × 070	087 × 070
P.A.	−873	758	768	768
rms (Jy beam−1)	0.022	0.016	0.009	0.008
Integrated intensity (Jy km s−1)	109 ± 11	46 ± 5	18 ± 2	5.8 ± 0.6
Peak flux (Jy beam−1)	0.86	0.87	0.45	0.23

Note. ^aThe velocity resolution indicated by the channel width is smaller than the true velocity resolution as the data are Hanning smoothed along the spectral dimension (Lundgren et al. 2012).

Download table as:  ASCIITypeset image

In order to produce spectral images with such high fidelity, we developed a CLEAN mask based on the emission pattern expected from a disk in Keplerian rotation (Beckwith & Sargent 1993), which requires some assumption about the stellar mass, disk inclination, disk size, phase offset, position angle (P.A.), and systemic velocity (see Section 4.1). Except for the stellar mass, all of these parameters can be estimated from an initial imaging of the data with a square mask; the final imaging results are consistent for any reasonable choices of these parameters. We identified the region of the rotating disk where the projected line of sight velocities corresponded to each individual spectral channel. This region was then convolved with a fixed Gaussian kernel (σ = 07) to broaden the mask and ensure that all of the significant emission was covered and that the CLEANed area was always much larger than the synthesized beam size. This recipe produces a separate mask for each channel, within which the measurement set was CLEANed deeply. Some additional CLEAN iterations were then applied with no mask. This masking technique significantly improves the quality of the images compared to the reference maps provided in the SV data package (the signal-to-noise ratio per channel improves by a factor of ∼2–3).

The ¹²CO J = 3–2 emission observed toward the disk around HD 163296 is asymmetric in its spatial morphology across the major axis of the disk (P.A. = 312°). For channels offset from the line center by velocities ≳ 1 km s⁻¹, the resolved emission appears systematically closer to the southern semi-minor axis (see Figure 2 and the top set of channel maps in Figure 1). Furthermore, channel-by-channel, the southern half of the disk is brighter than the northern. The ¹²CO J = 2–1 line has these same features (refer to Figure 1), but the effect is much less obvious; the morphological asymmetry can only be clearly seen for a few channels, |Δv| ≈ 1.7–2.0 km s⁻¹. The ¹³CO and C¹⁸O spectral emission have no apparent spatial asymmetries (see Figure 1).

Before attempting to explain the detailed morphologies of these CO data, we provide a brief primer on molecular line emission from a rotating disk (e.g., Omodaka et al. 1992; Beckwith & Sargent 1993). The observed line is centered around some transition frequency, f₀, with a linewidth that is sensitive to the quadrature sum of both the thermal velocity and turbulent motions of the gas. If the emitting gas has some bulk velocity along the line of sight, v_los, then the line center will be Doppler shifted away from the transition frequency, f₀, by an amount Δf = −v_los(f₀/c). Observations taken by (sub-)millimeter interferometers like ALMA have both spectral and spatial resolution, and produce an image of the source for a set of frequency, or equivalently velocity, channels. These spectral line observations comprise a three-dimensional data cube of position, position, and velocity. For an inclined disk in Keplerian rotation, the emission observed in a single channel "highlights" what part of the disk has that same projected velocity as it is Doppler shifted into the channel center.

The morphology of the line emission in any given channel changes depending on both the physical conditions and bulk motions of the gas. Modeling the line emission requires a calculation of the molecular excitation state of the gas and an integration of the radiative transfer equation along each line of sight, s, through the disk structure,

$$\begin{equation} I_\nu = \int _0^\infty S_\nu (s) \exp \left[-\tau _\nu (s)\right] K_\nu (s) ds, \end{equation} \tag{ 1 }$$

where $\tau _\nu (s) = \int _0^s K_\nu (s^{\prime })ds^{\prime }$ is the optical depth, K_ν is the absorption coefficient, and S_ν is the source function. All three of these terms depend upon the local temperature, T_gas, and density, ρ_gas, of the disk. Determining the structure of a disk using optically thick CO lines is difficult, and so previous studies have leveraged observations of multiple lines to break the model degeneracies (Dartois et al. 2003; Piétu et al. 2007).

To parse what these data reveal about the disk hosted by HD 163296, we will first present a toy model that qualitatively explains the observed morphological asymmetries noted in the ¹²CO transitions (Section 4.1). Next, we describe a more complex modeling procedure for testing whether a given disk structure produces the observed line intensity in addition to the morphological and brightness asymmetries (Section 4.2). We then present a simple model developed for pedagogical purposes (Section 4.2.1) followed by a more realistic hydrostatic disk model (Section 4.2.2). Finally, we demonstrate that some care must be taken when defining the rotation pattern of the disk, in light of the disk structure we derive with its large radial extent and molecular layer extending relatively high above the midplane (Section 4.3).

4.1. A Toy Model

The asymmetric shape of the ¹²CO J = 3–2 line emission cannot be explained by the emission pattern predicted for a vertically thin rotating disk. The coordinates of an inclined thin disk (x, y, z = 0), are related to the sky-plane relative to the observer, (x', y'), by

$$\begin{equation} \left(\begin{array}{@{}c@{}}x \\ y \\ z\end{array}\right) = \left(\begin{array}{@{}c@{}}x^{\prime } \\ y^{\prime } / \cos i \\ 0\end{array}\right), \end{equation} \tag{ 2 }$$

where i is the disk inclination (0° is seen face-on). For a disk rotating in circular, Keplerian motion, the projected line of sight velocity as a function of position is

$$\begin{equation} v_{\rm obs}(x^{\prime },y^{\prime }) = \sqrt{\frac{{GM}_\ast }{r}} \sin i \cos \theta, \end{equation} \tag{ 3 }$$

with an argument $\theta = \arctan (y/x)$ and modulus $r=\sqrt{x^2+y^2}$. This function is mirror symmetric across the x'-axis and so the shape of the emission observed at some velocity, v_obs, should appear symmetric about the x'-axis of the sky-plane. Of course, this geometrical argument is grossly simplified and considers neither the physical structure of the disk, nor the associated radiative transfer effects for spectral line emission (see Section 4.2). However, protoplanetary disks are thought to be geometrically thin (Pringle 1981), with radial surface densities and temperatures that roughly follow power-laws (i.e., slowly changing ρ and T in the outer disk; Weidenschilling 1977a; Adams et al. 1988). Under these conditions, this approximation works well and has been both a useful and effective guide in the study of line emission from protoplanetary disks. However, the observation analyzed here reveals an emission pattern that requires a modified conceptual framework (see also Pavlyuchenkov et al. 2007; Semenov et al. 2008; Guilloteau et al. 2012).

The toy model that we propose is a differentially rotating double cone, oriented along the z axis with an opening angle ψ ∈ (0, π/2) measured from the xy plane (i.e., the disk midplane). In this case, rays originating from the sky-plane will intersect the surface at

$$\begin{equation} \left(\begin{array}{@{}c@{}}x \\ y \\ z\end{array} \right) = \left(\begin{array}{@{}c@{}}x^{\prime } \\ \ y^{\prime } / \cos i + t \sin i \\ t\cos i\end{array}\right), \end{equation} \tag{ 4 }$$

where t solves the quadratic equation⁵

$$\begin{eqnarray} 0 &=& t^2 \left[\cos (2i) + \cos (2\psi)\right] \nonumber\\ && - 2 \sin ^2(\psi)[{x^{\prime }}^2 + {y^{\prime }}^2 \sec ^2(i) + 2 t y^{\prime } \tan (i)]. \end{eqnarray} \tag{ 5 }$$

The positive and negative roots of this equation correspond to the near and far halves of the double cone, respectively (where we define the near half as the one closer to the observer). Any given line of sight will intersect the near cone at larger y than for the flat disk, by an additional factor tsin i (assuming i < 180°). Therefore, if this near cone is in circular, Keplerian rotation, ${v_\theta }^2 ={GM}/\sqrt{x^2 + y^2}$ with no v_z component, the observed isovelocity contours will be systematically shifted to smaller y. The same argument follows for the far cone, except that isovelocity contours are systematically shifted to larger y. If i > 180° (which is the case for the HD 163296 disk; see Section 4.2), then the effects are the same except the solutions switch to describe the opposite half of the double cone. Larger values of ψ will amplify the differences between these curves. The emission morphology predicted from this double cone structure clearly mimics the observed asymmetry in the ¹²CO J = 3–2 line, as demonstrated in Figure 3.

Zoom In Zoom Out Reset image size

The brightness and morphological asymmetries of the ¹²CO J = 3–2 and ¹²CO J = 2–1 lines can be understood if the τ = 1 emitting surface is something like a double cone structure with ψ ∼ 15°. The bright portion of the line corresponds to the near cone that is systematically shifted along the minor axis (to smaller y). The dim "shadows," which are asymmetric in the opposite sense, represent the far cone. Figure 4 shows the velocity pattern of this toy model overlaid on the regridded channel maps of the ¹²CO J = 3–2 emission. The key feature of the double-cone is that it mimics an emission structure that has little contribution from the disk midplane: the structure is in some observational sense "hollow" at its center. Furthermore, the angle that the cone makes with the disk midplane (ψ ∼ 15°) is fairly large (but not abnormally so; van Zadelhoff et al. 2001) and suggests that the disk cannot be treated as geometrically thin.

Zoom In Zoom Out Reset image size

Detailed radiative transfer and chemical models (Aikawa & Herbst 1999; Aikawa & Nomura 2006; Semenov & Wiebe 2011; Walsh et al. 2012) predict a cone-like structure for molecules like CO. In these calculations, the midplane temperatures are coupled to the cold dust, greatly reducing the CO gas-phase abundance via freezeout onto dust grains. Higher in the disk, direct stellar radiation heats the material, liberating the CO from the grain surfaces and ensuring that CO is both abundant and warm. Above this molecular layer, where gas column densities are low, the CO is depleted by photodissociation via energetic photons from the star and background radiation field. It is this first layer, the cold midplane, that naturally explains the ALMA observations. In this region, the gas is too cold to densely populate the higher J CO transitions, so the observed emission is low. Additionally, the CO number densities should be heavily depleted by freezeout onto dust grains, further reducing the emission. Instead, most of the emission is from the warm molecular layer which, being high above the midplane, should produce asymmetric emission with respect to the major axis of the disk. This scenario also explains why the less abundant ¹³CO and C¹⁸O emission is not asymmetric: these lines originate from deeper layers in the disk, which correspond to lower ψ and less severe asymmetries that are not resolved by the observations. Indeed, calculations by Pavlyuchenkov et al. (2007) and Semenov et al. (2008) of HCO⁺J = 4–3 emission for axisymmetric, inclined disks (i ∼ 60°) predicted this same morphology.

4.2. Physical Model Analysis

While the toy model of a differentially rotating double-cone provides a useful qualitative explanation for the asymmetric morphology of the ¹²CO J = 3–2 line emission, we now undertake a more physically motivated, quantitative analysis. We build upon the disk structures developed by Dartois et al. (2003) and Qi et al. (2011) to reproduce the morphology of the ¹²CO J = 3–2 line along with the line intensities of the J = 2–1 emission from the ¹²CO, ¹³CO, and C¹⁸O isotopes. In this section, we introduce the modeling scheme that we will use to illustrate the sensitivity of the observations to the temperature structure and velocity field of the disk.

Following the formalism explained by Andrews et al. (2012), axisymmetric density and temperature structures are defined using a polar cylindrical coordinate system (r, z). The gas temperature structure, T_gas(r, z), is specified parametrically throughout this work. The gas surface density profile is assumed to be the self-similar model of a thin, viscous accretion disk (Lynden-Bell & Pringle 1974; Hartmann et al. 1998),

$$\begin{equation} \Sigma _{\rm gas}(r) = \Sigma _c \left(\frac{r}{r_c}\right)^{-\gamma } \exp \left[-\left(\frac{r}{r_c}\right)^{2-\gamma }\right], \end{equation} \tag{ 6 }$$

where r_c sets the radial size of the gas disk, γ is a power-law index, and $\Sigma _c = M_{\rm gas} (2-\gamma)/(2\pi r_c^2)$ is the normalization, where M_gas is the total gas mass. The density structure, ρ(r, z), is either parametrically defined (Section 4.2.1) or calculated from the temperature structure using the equation of hydrostatic equilibrium (Section 4.2.2). To calculate the CO number density, we assume that 80% of the gas (by number; μ = 2.37) is composed of molecular hydrogen and that the CO:H₂ relative abundance is a constant, f_co. Molecular depletion due to freezeout is implemented by decreasing f_co by a factor of 10⁸ in the region of the disk where gas temperatures fall below a freezeout temperature, T_gas < T_frz. A photodissociation boundary is calculated by vertically integrating the H nuclei density (0.706n_gas, the same procedure as in Qi et al. 2011; Aikawa & Herbst 1999) and thresholding for heights, z_phot, where the column density is less than a constant value,

$$\begin{equation} \sigma _s > 0.706 \int _{z_{\rm phot}}^\infty n_{\rm gas}(r,z^{\prime }) dz^{\prime }. \end{equation} \tag{ 7 }$$

This surface density is equivalent to the unitless Σ₂₁ defined by Aikawa & Nomura (2006) and used by Qi et al. (2011) where Σ₂₁ = σ_s/1.59 × 10²¹ cm⁻².

We assume a constant turbulent velocity width, ξ = 0.01 km s⁻¹, and disk inclination, i = 224°. The bulk gas velocities are described in each section and we fix several important model parameters: the stellar mass (M_* = 2.3 M_☉; Natta et al. 2004), source distance (d = 122 pc; van den Ancker et al. 1998), and major axis P.A. (measured east of north, P.A. = 312°). Note that the asymmetries in the ¹²CO J = 3–2 line emission effectively resolve the absolute disk spin orientation (Piétu et al. 2007), so our P.A. and inclinations differ by ∼180° from the models presented by Isella et al. (2007) and Qi et al. (2011). See Section 5 for a more detailed discussion.

For any given disk structure model, we calculate the CO level populations using the non-local thermodynamic equilibrium (non-LTE) molecular excitation and radiative transfer package LIME (Brinch & Hogerheijde 2010) and molecular data from the LAMDA database (Schöier et al. 2005; Yang et al. 2010). LIME is then used to generate channel maps matching the native channel spacing of the ALMA data (see Table 1). We then smooth using the Hanning algorithm in the spectral dimension with a three channel kernel to mimic the ALMA pipeline (see Section 5; Lundgren et al. 2012) and calculate model visibilities to match the Fourier sampling of the measured visibilities. We evaluate the model fit through a χ² statistic of the complex visibilities and also by inspection of CLEANed channel maps using the same masks and imaging procedure applied to the data (see Section 2).

It is important to emphasize that the models developed here involve many degenerate parameters: the exact structures we present are neither unique nor necessarily optimal. The analysis presented here is only meant to provide a basic physical understanding of some new, key aspects of these data.

4.2.1. A Pedagogical Structure

In this section, we will demonstrate how a vertical temperature gradient naturally produces the asymmetric emission morphology of the ¹²CO J = 3–2 line. In order to disambiguate the effects of density and temperature, we define a single parametric disk density structure and calculate the line emission for two simple cases for the temperature structure: (1) the disk is vertically isothermal and (2) the disk has a warm (T ∼ 30 K) molecular layer above a cold (T ∼ 20 K) midplane.

We assume that the gas is vertically distributed with the Gaussian profile appropriate for a vertically isothermal disk in hydrostatic equilibrium (neglecting self-gravity):

$$\begin{equation} \rho _{\rm gas}(r,z) = \frac{\Sigma _{\rm gas}(r)}{\sqrt{2\pi }H} \exp \left[-\frac{z^2}{2H^2}\right]. \end{equation} \tag{ 8 }$$

For this pedagogical example, the temperature and density structure are entirely decoupled. The disk scale height is defined to be H(r) = 16(r/150 AU)^1.35. The velocity field is defined such that the disk is in circular Keplerian rotation about the central star, with an azimuthal component

$$\begin{equation} {v_K}^2 = \frac{{GM}_\ast }{r}. \end{equation} \tag{ 9 }$$

Depletion due to CO freezeout is not considered in this calculation, but photodissociation is included with a threshold density σ_s = 5 × 10²⁰ cm⁻². The model disk has a gas mass of 0.09 M_☉, r_c = 115 AU, γ = 0.8, and f_co = 5 × 10⁻⁵.

We first calculate the ¹²CO J = 3–2 line emission from this disk structure assuming that it is vertically isothermal with a radial temperature structure set by a power law, T_gas = 65K (r/100 AU)^−0.5. The disk density and temperature structure is shown in Figure 5(a), and the corresponding channel maps in Figure 6 (top panels). The qualitative discussion from Section 4.1 is confirmed: the morphology of the line is essentially symmetric across the major axis with only a slight brightness asymmetry between the northern (near) and southern (far) sides of the disk (which is reversed if the inclination is flipped). This geometric effect, for sensitive data, indicates the absolute orientation of the disk and can qualitatively be understood by considering the emission that originates above the disk midplane. The half of the disk (located in the N) that is tilted toward the observer will then have a smaller projection on the observer sky plane (and appear dimmer) than the disk half (located in the S) that is tilted away. The double cone structure described in Section 4.1 similarly encapsulated this effect: in Figure 4 the isovelocity contours of the near cone (black curve) appear noticeably shorter in the N part of the disk than in the S for channels near the line center, |v_los| ≲ 1.02 km s⁻¹. Optical depth effects and the disk structure also contribute to the details of this brightness asymmetry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we build on that calculation by including a vertical gradient to the temperature structure. We start by defining a "warm" region where the gas is a constant T_a = 30 K for z > 20 AU (r/100 AU) or r < 150 AU (the freezeout radius for this disk; Qi et al. 2011). Outside this region (i.e., the disk midplane for r > 150 AU), the gas is a cooler T_m = 20 K. Since the gas high above the midplane and close to the star should be hotter than the 30 K typical for the outer region of this disk, we add a hot atmosphere above z_a = 73 AU (r/200 AU) with a power-law component T_a(r, z > z_a) = T_gas +40 K (r/200 AU)^−0.8. This disk structure is shown in Figure 5(b), and the corresponding channel maps in Figure 6 (bottom panels). Qualitatively, this model does an excellent job of producing both the morphological and brightness asymmetries noted in Section 3.

The disk structures described above are both massive and geometrically thick: the model scale height H is much larger than what the actual midplane temperatures would produce. Specifically, the defined H(r) would correspond to a midplane temperature of T_m = 50 K (r/100 AU)^−0.3, which is too warm to reproduce the line asymmetries (too much emission would be generated, essentially filling in the hollow cone emission morphology). However, the artificially large scale height is necessary to reproduce the morphology, since increasing the height of the emitting layer increases the projected separation of the front and back half of the disk (see Section 4.1). This suggests that a warm layer suspended above a cold midplane may help to raise the height of the emitting CO layer. It should be emphasized that a vertical temperature gradient is necessary to reproduce the data although it may not be the only contributing effect. Depletion of the gas phase CO due to freezeout onto grains should also reduce the contribution of the midplane emission, but this effect requires the cold temperatures that, by themselves, can produce this distinctive morphology.

4.2.2. A Hydrostatic Model Structure

While the previous model serves to demonstrate how a vertical temperature profile can explain the observed asymmetries in the ¹²CO J = 3–2 and ¹²CO J = 2–1 lines, the model structure is purely pedagogical and not physically self-consistent. We now present a hydrostatic disk model where the densities are internally consistent with the parametrically defined temperatures. This model roughly reproduces both the morphology and intensity of the observed emission simultaneously for the lines of interest, ¹²CO J = 3–2, ¹²CO J = 2–1, ¹³CO J = 2–1, and C¹⁸O J = 2–1.

We define the gas temperatures using a slightly modified version of the parameterization introduced by Dartois et al. (2003), where

$$\begin{equation} T_{\rm gas}(r, z) = \left\lbrace \begin{array}{@{}ll@{}}T_a + (T_m - T_a) \left[\sin {\frac{\pi z}{2 z_q}}\right]^{2\delta (r)} & \mbox{if z < z_q} \\ T_a & {\rm if}\ z \geq z_q \end{array}. \right. \end{equation} \tag{ 10 }$$

Here, T_m(r, z = 0) = 19 K (r/155 AU)^−0.3 is the temperature at the disk midplane and the atmosphere temperature above the height z_q is

$$\begin{equation} T_a(r, z) = 55\; {\rm K\,}(\sqrt{r^2 + z^2}/200{\rm \,AU})^{-0.5}. \end{equation} \tag{ 11 }$$

This prescription produces a smooth, monotonically increasing (with z) transition between T_m(r, z = 0) and T_a(r, z = z_q), with the vertical shape controlled by the parameter δ = 0.0034*(r − 200 AU) + 2.5 (with the requirement that δ ⩾ 0.3). The height of the disk atmosphere is assumed to have a radial distribution described by a truncated power law, z_q = 63 AU (r/200 AU)^1.3exp [ − (r/800 AU)²]. The midplane temperature profile we assume has the same CO freezout radius calculated by Qi et al. (2011), but with a shallower power-law index. The form of the vertical temperature profile (Equation (10)) was designed to mimic the dust temperatures in the disk models developed by D'Alessio et al. (1998, 1999, 2001, 2006). Since we do not do the radiative transfer calculation for a dust disk (e.g., Qi et al. 2011; Andrews et al. 2012), we adopt this parameterization as a natural choice for the outer regions of a disk.

The vertical distribution of the gas is determined by solving the differential equation for hydrostatic equilibrium with our parametrically defined temperature structure:

$$\begin{equation} -\frac{\partial \ln \rho _{\rm gas}}{\partial z} =\frac{\partial \ln T_{\rm gas}}{\partial z} + \frac{1}{{c_s}^2}\left[\frac{G M_\ast z}{(r^2 + z^2)^{3/2}}\right], \end{equation} \tag{ 12 }$$

where c_s² = k_BT_gas/μm_h is the sound speed. We normalize the density so that the total mass of the gas disk is M_gas = 0.09 M_☉, and assume γ = 0.8 and r_c = 150 AU. As before, we assume that the gas follows circular Keplerian orbits and that the self-gravity of the disk is negligible. We now include freezeout for T_frz = 19 K and photodissociation with σ_s = 9 × 10²⁰ cm⁻². For the CO isotopes, we assume the relative abundances measured for the interstellar medium by Wilson (1999): ¹²CO/¹³CO = 69 ± 6 and ¹²CO/C¹⁸O = 557 ± 30. The disk density and temperature structure are shown in Figure 7.

Zoom In Zoom Out Reset image size

In Figure 8 we present channel maps of predicted model ¹²CO J = 3–2 and C¹⁸O J = 2–1 emission, with the imaged residuals overlaid in color. This model successfully reproduces the distinctive spatial morphology of the ¹²CO J = 3–2 line (top panels) while also appearing symmetric in the C¹⁸O J = 2–1 line (bottom panels). The majority of the C¹⁸O emission is contributed from material much closer to the midplane than in the case of ¹²CO. The resulting signature of the temperature asymmetry is much less pronounced in the isotopologue lines, and is unresolved by these data. We show a comparison of the integrated line profiles for the data and model predictions in Figure 9 for the four CO lines of interest.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While this model matches the data fairly well, there are systematic residuals at large radii near the major axis at projected velocities of ≈1 km s⁻¹ where the model does not produce any significant emission (I_ν < 3σ). We cannot account for this emission by increasing the disk size (r_c), since emission at larger radii for channels at higher velocities as well as at the systemic velocity is overproduced. One interpretation of these significant residuals is that the gas at large radii is predicted to be moving too fast and so the emission morphology of the model appears to "fold" toward the major axis too slowly (when starting at the line center and looking channel by channel toward the blue or red shifted line wings). We address how the model may not accurately describe the true disk velocity field in the next section.

4.3. A Closer Look at the Gas Velocities

Although the models we present in Sections 4.2.1 and 4.2.2 leave non-negligible residuals compared to the data, they nevertheless serve to illustrate how a disk with a realistic vertical temperature structure naturally explains the asymmetric morphology of the ¹²CO J = 3–2 line. These data are an example of the exquisite sensitivity that ALMA will routinely achieve for observations of protoplanetary disks. With this in mind, we now explore how these data are sensitive to the bulk motions of the gas in the HD 163296 disk, necessitating the use of self-consistent, physically motivated disk structures in any detailed analysis of its structure or kinematics.

We begin with a broad discussion of the basic expectations for gas velocities in protoplanetary disks. These disks are usually assumed to be rotating differentially and in vertical hydrostatic equilibrium, so that there is no bulk vertical motion of the gas. At any location in the disk, there is an orbital velocity,

$$\begin{equation} \boldsymbol {v} = v \,\hat{\theta }, \end{equation} \tag{ 13 }$$

where $\hat{\theta }=(\sin \theta, -\cos \theta,0)$, and the gas is assumed to be on circular orbits. The standard (or tested, e.g., Dutrey et al. 1994) assumption is that v is equivalent to the Keplerian orbital velocity set by the gravitational potential of the central star,

$$\begin{equation} \frac{{v_K}^2}{r} = \frac{{GM}_\ast }{r^2}. \end{equation} \tag{ 14 }$$

Equation (14) implicitly assumes that the disk is geometrically thin, so that for a given r the gas at the midplane (z = 0) is moving at the same velocity as the gas suspended above it. But including the vertical geometry actually decreases the force from the star as well as the radial projection of that force vector (which sets the orbital velocity). This effectively reduces v as z increases, and implies that the disk has differential rotation in both the radial and vertical dimensions,

$$\begin{equation} \frac{v^2}{r} = \frac{r}{(r^2 + z^2)^{1/2}} \left(\frac{{GM}_\ast }{r^2 + z^2}\right). \end{equation} \tag{ 15 }$$

Recasting in terms of the toy model described in Section 4.1, Equation (15) reduces the observed velocity of the CO emitting layer by a factor of (cos ψ)^3/2 ∼ 0.95–0.91 for a cone angle, ψ, of 15°–20°. This rough estimate alone suggests caution when adopting a Keplerian, thin disk velocity field (Equation (14)) for modeling emission from molecules located high above the midplane. In principle, this assumption can be checked a posteriori for any disk model by calculating ψ (or equivalently z) of the τ = 1 emitting surface (see Section 5).

In addition to the gravitational potential of the central star, any radial change in the gas pressure will provide a force term and alter the gas velocities from the fiducial Keplerian field, v_K. Including the pressure gradient, the force equation becomes

$$\begin{equation} \frac{v^2}{r} = \frac{{GM}_\ast r}{(r^2 + z^2)^{3/2}} + \frac{1}{\rho _{\rm gas}}\frac{\partial P_{\rm gas}}{\partial r}. \end{equation} \tag{ 16 }$$

As both the gas temperature and density tend to decrease with r, a corresponding negative pressure gradient causes the gas to slow down and orbit at sub-Keplerian velocities. This phenomenon is thought to play an important role in the radial migration and growth of solids in these disks (Weidenschilling 1977b; Takeuchi & Lin 2002; Birnstiel et al. 2010), but is a subtle effect that slows down the gas from the fiducial Keplerian velocities by a rough factor of $(1-c_s^2/v_K^2)^{1/2} \sim 0.99$ for a temperature of 30 K at r ∼ 500 AU in the disk (Armitage 2009). However, if the density is falling faster with radius than a power-law, the pressure gradient will be larger and the gas velocities markedly slower. Therefore, the significance of this force term depends intimately upon the density and temperature structure. The popular self-similar solution for Σ_gas(r) that we have utilized features an exponential tail that will enhance this effect.

The last force term we consider is the self-gravity of the disk. Unlike the two previous terms which tended to slow down the gas and produce sub-Keplerian velocities, the additional mass contributed by the disk should increase v. The composite force equation is then

$$\begin{equation} \frac{v^2}{r} = \frac{{GM}_\ast r}{(r^2 + z^2)^{3/2}} + \frac{1}{\rho _{\rm gas}}\frac{\partial P_{\rm gas}}{\partial r} + \frac{\partial \phi _{\rm gas}}{\partial r}, \end{equation} \tag{ 17 }$$

where ϕ_gas is the potential due to self-gravity of the disk. The right-hand side of this equation includes three force terms: F_{stellar gravity}, F_{pressure gradient}, and F_{disk gravity} (appearing from left to right in Equation (17)). To self-consistently construct a disk structure now requires iteratively solving the equation of hydrostatic equilibrium,

$$\begin{equation} -\frac{\partial \ln \rho _{\rm gas}}{\partial z} = \frac{\partial \ln T_{\rm gas}}{\partial z} + \frac{1}{{c_s}^2}\left[\frac{{GM}_\ast z}{(r^2 + z^2)^{3/2}} + \frac{\partial \phi _{\rm gas}}{\partial z}\right], \end{equation} \tag{ 18 }$$

and Poisson equation,

$$\begin{equation} \nabla ^2 \phi _{\rm gas} = 4 \pi G \rho _{\rm gas}. \end{equation} \tag{ 19 }$$

For most disk configurations, the Poisson equation must be calculated numerically which we do using a fixed grid (Swarztrauber & Sweet 1975) without any iterations with the hydrostatic equilibrium equation (i.e., we do not include self-gravity when calculating the disk density structure for reasons explained below). As an illustration, we can roughly estimate the contribution of this correction by considering a special case with an analytic solution to the Poisson equation: a thin disk with a surface density profile Σ_gas = Σ₀r₀/r and no outer edge has

$$\begin{equation} \left.\frac{\partial \phi _{\rm gas}}{\partial r}\right|_{r=R} = 2\pi G\Sigma (R) \end{equation} \tag{ 20 }$$

at the midplane (Mestel 1963; Lodato 2007). In that scenario, the gas velocities are increased by a factor of $(1+ 2\pi Gr\Sigma (r)/ v_{\rm Kep}^2)^{1/2} \sim 1.005$–1.001 for radii ∼100–500 AU in the disk. We note that throughout our analysis we have only considered circular orbits. Including eccentricity in the disk (e.g., Regály et al. 2011) introduces a host of new parameters, asymmetries in the integrated line profile (which we have not investigated), and is beyond the scope of our analysis.

We now explore whether these ALMA data can distinguish between models that include these subtle effects on the gas velocities: (1) the differential rotation due to the vertical geometry of the disk, (2) the sub-Keplerian velocities of the gas due to the radial pressure gradient, and (3) the super-Keplerian velocities introduced by self-gravity. In order to disambiguate these three effects, we will consider separately the three force terms in Equation (17) that determine v: F_{stellar gravity}, F_{pressure gradient}, and F_{disk gravity}. When modeling the last two disk-specific terms, we revert to the thin-disk approximation where F_{stellar gravity} = GM_*/r². We use the same observing geometry, disk density, and temperature structure introduced in Section 4.2.2, changing only the gas velocities. As mentioned previously, the self-gravity of the disk affects both the velocity and density structures of the disk. Since both of these structures will impact the observed line emission and we are interested here in the effect of the former only, we choose to omit the disk potential in the equation of hydrostatic equilibrium and leave the disk density structure unaltered. For completeness, we also calculate the emission from a model whose velocities are determined by all three terms (Equation (17)). Due to the more complicated velocity field and the fine spatial sampling of our disk models, we found that the packaged raytracer in LIME was prohibitively slow. Therefore, for these models we only used LIME to calculate the non-LTE level populations. To generate the model images, we utilized the axisymmetry of our disk models to interpolate onto a fine two-dimensional grid (Fan et al. 2005) and integrated the radiative transfer equation. None of these four models require any additional parameters.

Figure 10 shows the fractional and absolute velocity difference compared to the fiducial Keplerian velocity field (v_K² = GM_*/r) for each model structure. As expected, the geometry of the disk height and the disk pressure gradient slow down the gas high above the midplane and at large disk radii (r ≳ 300 AU). The disk self-gravity speeds up the gas, but has a smaller effect than the first two terms. Combining all three effects for our fiducial model, the fractional difference in the gas velocities are on the order of a few percent in the midplane. However, the absolute velocity difference in the CO emitting region in the outer disk can be as great as ∼0.1 km s⁻¹, the same as the native channel width of the Band 7 observations.

Zoom In Zoom Out Reset image size

The impact of these changes to the velocity field are significant at the ≳ 3σ level across many channels. Figure 11 shows the difference between the model with the standard Keplerian rotation (v² = GM_*/r; Section 4.2.2) and the individual models whose velocity fields are summarized in Figure 10. For all of these models, the central channels were the least affected since the projected velocity (for this inclination) is small (v_los∝cos θ). The model that accounts for the vertical geometry of the disk differs most significantly, particularly for the higher velocity channels (v ≳ 1.4 km s⁻¹). The radial pressure gradient impacts the emission at large radii, where the surface density profile falls off sharply. By itself, the disk self-gravity has the least significant effect. However, the combined model is clearly not dominated by any one term, indicating that models for data of this high quality need to incorporate all of these more physical and self-consistent calculations for the velocity field. The usual assumptions can, of course, be instead tested a posteriori as we have done here.

Zoom In Zoom Out Reset image size

We have conducted an analysis of multiple CO emission lines observed toward the protoplanetary disk hosted by HD 163296. The exquisite spatial resolution and sensitivity of these ALMA SV observations, combined with good spectral resolution of the ¹²CO J = 3–2 line, clearly reveals that the spatial morphology of the emission is asymmetric across the major axis of the disk. We interpret this asymmetry as a resolved signature of the vertical temperature gradient in the disk generically predicted by physical/chemical models of disk structures (e.g., Pavlyuchenkov et al. 2007; Semenov et al. 2008). We developed a series of physically self-consistent disk structures, showing that these ALMA data can distinguish between models that correct the gas velocities for the vertical thickness of the disk and the radial pressure gradient of the gas.

We started with a toy model that can conceptually explain the asymmetric morphology of the ¹²CO J = 3–2 line. An emitting surface described by a rotating double cone that makes an angle ψ ∼ 15° with the disk midplane naturally mimicked the observed morphology of this disk. Such a toy model is analogous to a layered disk structure, where the contribution of the midplane emission is small due to the cold gas temperatures and depleted abundances produced by the freezeout of gas phase CO onto grain surfaces. We confirmed the qualitative behavior of this toy model by calculating the non-LTE emission from both pedagogical and hydrostatic axisymmetric disk models, concluding with a disk structure that roughly reproduced the observed morphology of the ¹²CO J = 3–2 emission. Our analysis also showed that a vertical temperature gradient, which has been reported in this disk by both Qi et al. (2011) and Akiyama et al. (2011), was consistent with the observed intensities of the J = 2–1 emission from ¹²CO, ¹³CO, and C¹⁸O, as well as the symmetric morphology of the rarer isotopologue lines.

Spectral line emission, even when it is spatially resolved, cannot easily define the absolute orientation of the disk: there is a 180° ambiguity in the inclination and P.A.. However, the morphological and brightness asymmetry highlighted here resolves the absolute orientation of the disk. We report values of the disk inclination and P.A. that differ by 180° from previous analyses based on less sensitive data (Isella et al. 2007; Hughes et al. 2008; Qi et al. 2011). For HD 163296, an independent check on the viewing geometry is provided by its outflow, where Doppler shift measurements of Hα, [S ii], [O iii], and [N ii] emission lines and imaging can differentiate between the jet (blue shifted, seen on the near side of the disk) and counter jet (red shifted, partially hidden behind the disk). For HD 163296, the blueshifted jet is oriented in the SW direction while the redshifted counter jet is moving NE (in agreement with the disk wind observed via ¹²CO by Klaassen et al. 2013). This is consistent with our reported inclination axis, which tilts the NE part of the disk toward the observer. Furthermore, the projection of the disk rotation axis on the observer sky plane (42°) is aligned with the P.A. of the jet (measured for the counter jet as 420 ± 35; Grady et al. 2000). The jet inclination ($51_{-9}^{+11} {}^\circ$; Wassell et al. 2006; Günther et al. 2013) also agrees with the disk inclination within the 1σ error bars. We confirm that the jet appears to be orthogonal to the disk plane, although the large error bars in the jet inclination leaves open the possibility that the jet is slightly skewed (Δi ≲ 10°).

A subtle, but important, point is that spatially resolving the asymmetry caused by the vertical temperature gradient requires that the data has both excellent spatial and spectral resolution. In Figure 12 we show a pixel map of the angle, ψ = tan ⁻¹(z/r), that the τ = 1 surface of our fiducial model (Section 4.2.2) makes with the disk midplane for the ¹²CO J = 3–2 and ¹²CO J = 2–1 emission. Both of these lines originate at nearly the same height, and both models exhibit an asymmetric emission morphology at a Δv = 0.1 km s⁻¹ resolution. However, when imaged at the ALMA SV resolutions, the J = 2–1 line appears mostly symmetric—just like the data (except for the slight rotation at |v| ≈ 1.7 km s⁻¹). Figure 13 shows the accumulated degradation of the spectral resolution for ¹²CO J = 2–1 model emission. The spectral averaging (middle panel), Hanning smoothing (right panel), and Fourier sampling (bottom row) all diminish the apparent morphological asymmetry even though the first two effects are operations along the spectral dimension. The muted asymmetry of the ¹²CO J = 2–1 emission can be attributed to the larger beam size and coarser spectral resolution of the Band 6 observations and not to the disk structure (as is the case for the rarer ¹³CO and C¹⁸O isotopologues). The Band 6 data have been averaged and smoothed in the spectral dimension (see Section 4.2), which smooths the correlated emission in the spatial dimension (Thompson et al. 1986). The consequence of a coarse spectral resolution is worth emphasizing, since in practice it may significantly degrade the utility of the excellent spatial resolution of ALMA.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explored the sensitivity of the Band 7 observations to the gas kinematics with a series of models where only the bulk velocities were changed. We found that these data could distinguish between models for which those gas velocity differences were ≳ 5%, corresponding to classes of models that incorporated the vertical geometry of the disk as well as the radial pressure gradient in the gas. Indeed, a few disks have been identified as having sub-Keplerian velocities (Wang et al. 2012), and we suggest that these phenomena provide an intuitive explanation (although might not be sufficient for such systems). For our fiducial model, the disk self-gravity did not appear to significantly alter the velocity field, but may be important when modeling the density structure. Furthermore, the contributions from the gas pressure gradient and self-gravity are model dependent, and so may become even more important for certain systems (e.g., Cesaroni et al. 2005; Bergin et al. 2013).

There are a few aspects of the disk hosted by HD 163296 that accentuate the subtle effects we have addressed. First, the disk is both large and relatively nearby so that the longest baselines in this SV dataset (∼400 kλ) correspond to spatial scales that will only be attained by full ALMA for T Tauri disks in more distant star-forming clusters. Second, HD 163296 is an A star that can passively heat its disk out to large radii. And so in addition to having a large radial extent and thick vertical geometry, the disk is warm and bright for sub-millimeter observations. Last, the intermediate disk inclination helps to evenly project both the vertical and radial dimension of the disk on the sky plane (Semenov et al. 2008). The inclination also projects a significant fraction of the intrinsic gas motions along the line of sight. This reduces the smoothing effect that the spectral resolution has upon the spatial pattern of the line (see Figure 13).

With high-resolution, sensitive instruments such as ALMA and the JVLA, there is a host of opportunities for studying circumstellar disks and phenomena associated with planet formation (e.g., Wolf & D'Angelo 2005; Semenov et al. 2008; Cossins et al. 2010; Cleeves et al. 2011; Gonzalez et al. 2012; Ruge et al. 2013). These observations demonstrate the necessity of more sophisticated, physically self-consistent approaches when analyzing data from this new generation of observing facilities.

We have analyzed sensitive, sub-arcsecond resolution observations of the ¹²CO J = 3–2, ¹²CO J = 2–1, ¹³CO J = 2–1, and C¹⁸O J = 2–1 emission lines from the protoplanetary disk hosted by HD 163296. The key conclusions of our analysis are the following.

• 1.  
  The ¹²CO J = 3–2 spectral line features a clear and systematic morphological and brightness asymmetry across the major axis of the disk. The ¹²CO J = 2–1 line exhibits similar, but muted, behavior. No asymmetries are observed for the rarer ¹³CO and C¹⁸O isotopologues.
• 2.  
  The resolved morphological and brightness asymmetries of the ¹²CO J = 3–2 line emission as well as the symmetric emission of the rarer isotopologues are the signatures of a vertical temperature gradient in the disk. A double cone structure encapsulates the salient physical features of a disk with a cold midplane and warm atmosphere and mimics the distinctive emission morphology by moving the surface that dominates the emission above (and below) the disk midplane.
• 3.  
  We presented a series of simple disk structures that demonstrated how the observed data can be roughly reproduced by a disk with a cold midplane and warm atmosphere. For our fiducial model, the τ = 1 surface for the ¹²CO emission was suspended ∼15° above the midplane, agreeing with our morphological characterization using the toy double cone structure.
• 4.  
  The excellent spatial and spectral resolution of the Band 7 observations are surprisingly sensitive to the bulk velocities of the gas and can distinguish between models that include the vertical geometry of the disk and radial pressure gradient (a fractional difference in the bulk gas velocity field of ≳ 5%). The inclusion of self-gravity for our fiducial disk model is a less important correction, but may become dominant for other disk models.
• 5.  
  The excellent sensitivity and spatial resolution of these ALMA SV data require careful processing of models. In particular, coarse spectral resolution of the data can strongly impact the spatial morphology of the observed emission.

We are grateful to Matt Payne, Diego Muñoz, Eugene Chiang, David Knezevic, Moritz Günther, Marc Metchnik, and Kees Dullemond for insightful conversations. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00010.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. We acknowledge support from NASA Origins of Solar Systems grant No. NNX11AK63.