ALMA OBSERVATIONS OF HD 141569's CIRCUMSTELLAR DISK

The dust emission is clearly resolved by the ALMA beam. The continuum (with the CO channels removed) is deconvolved and imaged using CASA's CLEAN algorithm. The average wavelength across the frequency range is 870 μm. A threshold of $\tfrac{1}{2}\times {\sigma }_{\mathrm{rms}}$ and a natural weighting are used to produce the final cleaned product in Figure 1 (with contours corresponding to 3, 6, 12, and 21 × σ_rms). The inner disk around HD 141569 is imaged out to 56 au (assuming a distance of 116 pc). The longest baseline is unable to resolve a central clearing of <15 au, leading to a central peak near the pointing center (the star and inner disk). The peak intensity in the cleaned data is 1.74 mJy beam⁻¹, corresponding to an S/N of ∼25. At 870 μm, the thermal emission from the host star's photosphere contributes <1% to the peak flux per beam, assuming a blackbody with T_Eff = 10,500 K, a radius of 1.7 R_⊙, and a distance of 116 pc. The star's flux is thus negligible, as long as corona and chromospheric effects can be ignored.

Zoom In Zoom Out Reset image size

The dust distribution is constrained using CASA's uvmodelfit, which fits single component models directly to the visibility data and selects the best fit through χ² minimization. We run this task to fit a uniform disk model to the continuum data (the CO channels are split out) and list the best-fit model in Table 3. Disks with inclinations near i ∼ 55° are favored with a major axis of about 085, corresponding to ∼85 au at a distance of 116 pc. The preferred model has a total continuum flux density of 3.78 ± 0.23 mJy. This is within 15% of the flux density found by summing the total flux from the cleaned image down to the 3σ contour. The uncertainty in the flux is dominated by the uncertainty in the absolute flux scale, which is taken to be 10%. This sets our flux estimate of the inner dust disk to be 3.8 ± 0.4 mJy.

Table 3.  Summary of CASA's uvmodelfit Results for the Debris Disk

Parameter	Continuum (Debris)
Flux Density	3.78 ± 0.23 mJy
X Offset	−0032
Y Offset	−0023
Major Axis	085
Axis Ratio (inclination)	0.58 [55°]
Position Angle	−88

Note. The data were fit by comparing a simple, uniform disk model to the data visibilities. The fitting uncertainties for parameters other than flux are not included here, but are addressed for the gaseous disk in Section 5.

Download table as:  ASCIITypeset image

In addition to the continuum, CO(3-2) emission is kinematically and spatially resolved using the FDM capabilities of the ALMA correlators, with a spectral resolution of 0.85 km s⁻¹. The double-horned spectrum is shown as a function of LSRK velocity in Figure 2. The previously constrained system velocity of 6 km s⁻¹ is shown, as well as the asymmetric emission from the disk (Dent et al. 2005).

Zoom In Zoom Out Reset image size

The CO is continuum subtracted using the CASA task uvcontsub. Figure 3 (left panel) shows the brightness map for the CO line (zeroth moment), in which the 3σ CO contour extends out to 18 (∼210 au). This is compared directly with the continuum emission (contours), which is more centrally concentrated. The right panel shows the velocity map (first moment), with the CO brightness contours overlaid. There are two brightness peaks, each at about ∼0.9 Jy km s⁻¹. The peaks are separated by ∼05 in a morphology that resembles ring ansae and is suggestive of an inner gas cavity. There is only a tenuous CO detection within this ∼05 (∼50 au) diameter cavity, which is broadly consistent with previous shorter wavelength observations that find only tenuous CO between about 10 and 50 au in radius (Brittain & Rettig 2002; Goto et al. 2006).

Zoom In Zoom Out Reset image size

The velocity field map shows clear Keplerian rotation, with the gas south of the star approaching us. The brightness is skewed westward (right in the image), relative to the velocity map, which is discussed in more detail below.

Figure 4 shows maps for 25 velocity channels between −0.5 and 11.5 km s⁻¹. Contours represent 3, 6, 9, and 24 times the rms noise of the zeroth moment. The spectral resolution of the velocity, 0.85 km s⁻¹, is a factor of two larger than the channel width. For Figure 4, velocity channel spacing is chosen to be 0.50 km s⁻¹ to include a slight oversampling. The total flux density of the CO given in Table 2 is determined by summing the flux in the zeroth moment map down to 3 × σ_rms and multiplying by the number of beams. This value is consistent with integrating over all channels of the CO map to within 10%. Note again that there is a clear asymmetry in the emission west of the star.

Zoom In Zoom Out Reset image size

The peak flux in the northwestern limb is significantly brighter than its counterpart in the northeastern and southwestern limbs.

As shown in Figure 3, the high spatial and velocity resolution capabilities of ALMA yield a well-constrained velocity field. These data can thus be compared to a Keplerian disk model to infer system properties. Trial models are generated by first assuming a uniform Keplerian disk. For simplicity, the inner cavity, temperature profile, and line broadening of CO (which is expected to be small) are not factored into the model. Each model is projected to the disk geometry and the LSRK velocity is subtracted. The model is then convolved with a 2D Gaussian beam as given in Table 2. Using MCMC techniques (specifically, Metropolis–Hastings with Gibbs sampling), the posterior distributions are calculated for the disk's inclination, position angle, LSRK system velocity, dynamical center, and mass. We assume flat prior distributions over the ranges given in Table 4. Model comparison is conducted in the image domain due to the high-velocity resolution and signal-to-noise.

Parameter space is explored through a random walk directed by Metropolis–Hastings MCMC (e.g., Ford 2005). For each new trial, two model parameters are randomly chosen and then updated by drawing a Gaussian random parameter centered on the current model (state i). The acceptance probability for the new trial model (state i + 1) is given by

$$\begin{eqnarray}&&\alpha ={\rm{\min }}({e}^{\displaystyle \frac{1}{2}({\chi }_{i}^{2}-{\chi }_{i+1}^{2})},1),\end{eqnarray} \tag{ 1 }$$

where we take

$$\begin{eqnarray}&&{\chi }_{i}^{2}=\sum \displaystyle \frac{{(D-{M}_{i})}^{2}}{{\sigma }^{2}}.\end{eqnarray} \tag{ 2 }$$

Here, D is the data from the CO first moment map (see Figure 2), M_i is the current model, and σ = 0.5 km s⁻¹ is the velocity channel width. The summation is over all points on the moment map. If α is greater than a random number drawn from a uniform [0,1] distribution, then the new model is accepted and recorded in the Markov chain. If the model is rejected, then the previous model is used again and re-recorded.

The MCMC routine is run using three chains, each with randomly chosen starting points in the flat prior parameter space. Each chain contains 100,000 links of which about 1000 are needed for burn-in. The three chains converge on similar parameters, and the distributions are combined to give the resulting posterior distributions in Figure 5. The blue points correspond to the values of highest probability. The most probable parameters (i.e., the mode of the distributions) are given in Table 5. Uncertainties are given by a 95% credible interval unless otherwise stated. The most probable mass is 2.39 M_⊙, for a distance⁹ of 116 pc. Since there is a degeneracy in inclination and mass, we give Msin(i) and M. Previously constrained stellar mass estimates are between 2.0 and 3.1 M_⊙ (e.g., Merín et al. 2004; Wyatt et al. 2007). The posterior distributions for both quantities are sampled independently by the MCMC. Ultimately, the uncertainty in the derived mass is dominated by the distance uncertainty. The re-analyzed Hipparcos Catalog distance with 1σ uncertainty is 116 ± 8 pc (Van Leeuwen 2007). Considering only this 1σ distance uncertainty with our most probable mass yields 2.39^+0.16_−0.16 M_⊙.

Zoom In Zoom Out Reset image size

The most probable parameters are used to construct a final disk model, which is shown in Figure 6. The residuals of the model are also shown as percent deviation from the data. The most probable model typically shows agreement with the data to about 10%, but has larger deviations along the minor axis of the data/model.

Zoom In Zoom Out Reset image size

An interesting asymmetry is observed in the CO channel map. Looking at the "butterfly" features in the 4–8 km s⁻¹ channels, there is a localized flux enhancement on the northwestern (top right) component of the gas. The east wing of the butterfly has a fairly symmetrical intensity about the system velocity of 6 km s⁻¹, while the west wing is asymmetrical.

To explore this feature further, Figure 7 shows three of the channel maps (4.5, 6, and 7.5 km s⁻¹), along with the continuum using 3, 6, 9, and 12 × σ_rms contours. The southern components of both sides appear to be approximately symmetric, but a strong asymmetry becomes obvious for the 6 and 7.5 km s⁻¹ maps, in which the western wing is brighter than the eastern wing by ∼40% in each channel. These channel maps also suggest that there is indeed an inner cavity to the CO disk, as noted in other studies (Goto et al. 2006; Flaherty et al. 2016).

Zoom In Zoom Out Reset image size

Since the asymmetry is present throughout multiple channels (see Figure 4), the feature appears to be real in the data. While the exact source of the flux enhancement is unknown, it may be caused by asymmetries in the inner disk edge, such as vortex formation (e.g., Lyra et al. 2008a, 2008b) or by perturbations from an unseen companion. Dent et al. (2014) observe a large asymmetry in β-Pic that is attributed to localized collisions of gas-rich comets. The asymmetry is also in the general direction of the two distant red dwarf companions that orbit at ∼1000 au. Follow-up observations and detailed simulations are required to determine the cause of the CO disk morphology.

An initial estimate for the dust mass is made by assuming that the emission is optically thin, dominated by millimeter grains, and spatially concentrated in a thin ring. In this case,

$$\begin{eqnarray}&&M=\displaystyle \frac{4}{3}{\rho }_{i}\pi {s}^{3}\displaystyle \frac{{F}_{\nu }(\mathrm{Obs})}{{B}_{\nu }(R){{\rm{\Omega }}}_{{\rm{s}}}},\end{eqnarray} \tag{ 3 }$$

where F_ν(Obs) is the observed flux density of the continuum, B_ν(R) is the blackbody intensity for a single grain placed at a distance R from the star, and Ω_s is the solid angle of a single grain. The grains are further assumed to be in thermal equilibrium with the host star, to have an internal density ρ_i = 2.5 g cm⁻³ and size s = 1 mm, and to be perfect absorbers and radiators (albedo of 0, emissivity of 1). We note that this mass estimate is equivalent to $M=\tfrac{{d}^{2}{F}_{\nu }(\mathrm{Obs})}{{\kappa }_{\nu }{B}_{\nu }(R)}$ with ${\kappa }_{\nu }=3\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, for our assumptions. This opacity is within a factor of two of the millimeter opacity used by Flaherty et al. (2016).

For an approximate lower limit, the ring can be envisaged to be at R = 10 au, which represents the innermost location for large grains based on SED modeling (Malfait et al. 1998). At this distance and for the noted assumptions, the grains would be T ∼ 200 K, which yields a millimeter grain mass¹⁰ of $0.04\,{M}_{\oplus }$. Placing grains at larger stellar separations would require additional mass to explain the emission. For example, if all the grains were placed at R = 50 au (T ∼ 90 K), the millimeter grain mass would be ∼0.09 M_⊕.

This simple estimate may only correspond to the actual dust mass if the observed millimeter grains are leftovers that were never incorporated into planets. Instead, if the grains are produced by the evolution of a nascent debris disk, the total mass can be significantly different. We explore this possibility next using a size distribution of grains spread throughout a disk.

For simplicity, we assume that the surface density of material decreases as r⁻¹ over disk radii 10–50 au. The dust is assumed to be absent outside of these boundaries. We further take the grains to radiate efficiently as long as the their diameter (2s) is equal to or larger than the absorbing/emitted photons (e.g., Wyatt & Dent 2002). For wavelengths larger than the grain's diameter, the emission and absorption coefficients (${Q}_{\nu }({\rm{em}})={Q}_{\nu }({\rm{abs}})={Q}_{\nu }$) are inversely proportional to the photon wavelength. Specifically,

$$\begin{eqnarray}{Q}_{\nu }=\left\{\begin{array}{lc}1 & \,2s\gt \lambda \\ \tfrac{2s}{\lambda } & \,\mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

We only consider a "total" debris mass up to some maximum parent body size, which is taken to be s_max = 50 km. This does not mean that 50 km is envisaged to be the largest solids in the debris disk; it is only the maximum size we consider in a given size distribution. To get a total debris mass for solids s < 50 km, a particle size distribution must be assumed. Lacking further constraints, we use a collisional cascade such that the mass per size increment ${M}_{s}\propto {s}^{-0.5}$ (e.g., Dohnanyi 1969). The total mass is then determined by requiring the model continuum flux density to match our observations. In practice, the debris disk is divided into a series of rings (here 100), placed evenly between 10 and 50 au. If each ring has the same mass, then the surface density profile follows ${r}^{-1}$. A flux density for each ring is then calculated by first deriving a grain temperature, assuming that the grains are dark (albedo ∼ 0) and balancing the received and emitted powers using a blackbody model with the effects of Q_ν. The grain temperature (e.g., Wyatt & Dent 2002) is

$$\begin{eqnarray}&&{T}_{{\rm{g}}}={T}_{{\rm{g}},\mathrm{BB}}{\left(\displaystyle \frac{{Q}_{{\rm{abs}}}({T}_{\mathrm{star}})}{{Q}_{{\rm{abs}}}({T}_{{\rm{g}}})}\right)}^{1/4},\end{eqnarray} \tag{ 5 }$$

where ${T}_{\mathrm{star}}$ is the host star's surface temperature (assuming it is a blackbody), ${T}_{{\rm{g}},\mathrm{BB}}$ is the equilibrium grain temperature if the grain were also a perfect blackbody, and T_g is the actual grain temperature. The equation must be solved iteratively, but converges quickly. For HD 141569, ${T}_{{\rm{g}}}\approx {T}_{{\rm{g}},\mathrm{BB}}$ except at grains less than 10 s of microns. To calculate Q_abs(T), Q_ν is integrated over all frequencies and weighted by a blackbody of the given temperature, i.e.,

$$\begin{eqnarray}&&{Q}_{\mathrm{abs}}(T)=\displaystyle \frac{\int {B}_{\nu }(T,\nu ){Q}_{\nu }d\nu }{\int {B}_{\nu }(T,\nu )d\nu }.\end{eqnarray} \tag{ 6 }$$

Taking T_star = 10,500 K and the above grain size and spatial distribution, we find that M(s < 50 km) ∼ 160 ${M}_{\oplus }\tfrac{{\rho }_{i}}{2.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}}$.

This result should be interpreted with caution. A steeper (shallower) solid size distribution can lead to significantly larger (smaller) masses. The result is also dependent on the internal density of the grains, as well as their effective albedo and emissivity. Nonetheless, the results are illustrative that significant debris may be distributed between 10 and 50 au. The total mass of solids would be much larger should debris (at a lower surface brightness) be present at disk radii $r\gt 50\,{\rm{au}}$, which would be consistent with single-dish measurements (see Table 1).

The mass of an optically thin gas disk near LTE can be calculated from the integrated line intensity (e.g., Perez et al. 2015). Given a line flux of F_OBS = 15.7 Jy km s⁻¹, the average line intensity over the source's solid angle Ω is

$$\begin{eqnarray}&&\hat{I}=\displaystyle \frac{{F}_{{\rm{OBS}}}}{\lambda {\rm{\Omega }}},\end{eqnarray} \tag{ 7 }$$

where λ = 867 μm is the average wavelength of the observations. The upper transition level column density of CO is given by

$$\begin{eqnarray}&&{N}_{3}=\displaystyle \frac{4\pi \hat{I}}{h\nu {A}_{32}},\end{eqnarray} \tag{ 8 }$$

where ν = 345.79 GHz is the frequency of the molecular feature, and A₃₂ = 2.497 × 10⁻⁶ Hz is the Einstein absorption coefficient¹¹ for the transition.

In the following, J = 3 (the upper transition level) unless otherwise noted (such as in the summation). Under the assumption that all J energy levels are populated in LTE, the total column density is given by

$$\begin{eqnarray}&&{N}_{{\rm{Total}}}={N}_{J}\displaystyle \frac{Z}{2J+1}{e}^{\displaystyle \frac{{{hB}}_{e}J(J+1)}{{kT}}},\end{eqnarray} \tag{ 9 }$$

and Z is

$$\begin{eqnarray}&&Z=\displaystyle \sum _{j=0}^{\infty }(2j+1){e}^{-\displaystyle \frac{{{hB}}_{e}j(j+1)}{{kT}}}.\end{eqnarray} \tag{ 10 }$$

Here, B_e = 57.635 s⁻¹ is the rotational constant⁴, T is the gas temperature, Z is the canonical partition function. The gas mass is then given by

$$\begin{eqnarray}{M}_{{\rm{CO}}} & = & {m}_{{\rm{CO}}}{N}_{{\rm{Total}}}{\rm{\Omega }}{d}^{2},\\ & = & \displaystyle \frac{4\pi {m}_{{\rm{CO}}}\,{d}^{2}\,{F}_{{\rm{OBS}}}\,Z}{h\nu \lambda {A}_{32}(2J+1)}\,{e}^{\displaystyle \frac{{{hB}}_{e}J(J+1)}{{kT}}}\end{eqnarray} \tag{ 11 }$$

for a solid angle Ω and distance to the object d. Taking a gas temperature of T = 33 K, the minimum excitation temperature of the J = 3–2 line, gives ${M}_{{\rm{CO}}}=1.9\pm 0.2\times {10}^{-3}\,{M}_{\bigoplus }$, with the uncertainty propagated from the CO flux density uncertainty in Table 2. The corresponding spatially averaged column density, N_Total, is 1.2 ± 0.1 × 10¹⁶ cm⁻². This is within the optically thin limit (Wyatt et al. 2015), but should not be taken as an independent confirmation because we assumed the gas to be thin for the mass calculation. If the gas is partly optically thick, then the actual CO gas mass could be larger by a factor of a few (Matrà et al. 2015). As such, the CO mass here could be interpreted as a lower limit. Due to the uncertainty in the appropriate amount of the gas, we will only report the CO mass as M_CO ∼ 2 × 10⁻³ M_⨁ to emphasize that the calculation has important unknowns.

Flaherty et al. (2016) find a gas model with total mass of 13⁺⁵⁰₋₉ M_⨁ as constrained by LTE models of gas temperature and density of CO(1-0) and CO(3-2) with CARMA and SMA, respectively. If we assume the 10⁴ ISM number density abundance ratio for H₂ to CO (as in Flaherty et al.), the inferred H₂ gas mass from the ALMA observations is ${M}_{{{\rm{H}}}_{2}}$ ∼ 1.4  M_⨁. Including additional metals would increase the total inferred gas mass to be slightly above ∼1.5  M_⨁, which is a factor of a few below the lower bound of the SMA and CARMA based model. The observations and models altogether thus suggest that there is one to a few tens M_⊕ of gas mass, assuming the ISM scaling can be used, which is not obviously the case. Additional caveats for these gas-mass estimates are discussed below.

The morphology of HD 141569 shows a dust disk extending out to about 56 au and an extended CO gas component between about 30 and 210 au. This structure alone suggests that the system is an evolved transition disk. However, as discussed below, HD 141569 may be better interpreted as a nascent debris system. The distinction is that the dust would be second generation, and any associated size distribution would reflect the clearing stages of planet formation rather than grain growth outcomes.

4.4.1. Primordial versus Second Generation