Modeling the Spatial Distribution and Origin of CO Gas in Debris Disks

"Lucy–Richardson deconvolution" is an iterative rectification method for observed distributions, presented independently by Richardson (1972) and Lucy (1974). This method attempts to restore the original distribution from an observed distribution, in which the latter corresponds to a degraded version of the former. The key to this method is that the original distribution can be recovered if sufficient a priori information about the degrading process is available. Examples of degrading processes are spatial distortion by an instrumental point-spread function (PSF), instrumental broadening of spectral lines (equivalent to a velocity PSF), additive noise, or other more complex processes (e.g., Zech 2013; Stock et al. 2015; Zorec et al. 2016).

The use of a priori kinematical information for deriving CO emissivity distribution in galaxies was demonstrated over 20 years ago by Scoville et al. (1983). By applying the Lucy (1974) iterative rectification scheme, Scoville et al. derived a de-projected emissivity distribution that successfully reproduced the observed line intensity profiles (in their case, the velocity field of the galaxy was known a priori via optical line studies). The advantage of this method is that it does not need any previous assumption on the surface brightness distribution (e.g., whether it is a power law or other functional form), and that it converges relatively rapidly (typically in less than ten iterations).

In the case of astrophysical disks (circumstellar or galactic), the observed line intensity profiles result from the (double) convolution between the spatial and velocity PSFs, with the intrinsic emissivity distribution. If x and y are the linear displacement coordinates parallel and perpendicular to the disk's major axis, we consider a point (x, y) located at distance R from the center of a circumstellar disk. If the velocity field is known at all points (e.g., Keplerian), ρ(R) is the axisymmetric emissivity distribution, and i the disk inclination angle, then in the case where there is no instrumental broadening the measured intensity in that pencil beam point is

$$\begin{eqnarray}&&J(x,y;v)=\displaystyle \frac{\rho (R)}{\cos i}\delta (v-{v}_{x,y}),\end{eqnarray} \tag{ 1 }$$

where v_{x, y} is the Keplerian velocity at position (x,y). In reality, the observed line intensity profiles correspond to convolution of the emissivity distribution ρ(R) with the instrumental PSF:

$$\begin{eqnarray}&&I(x,y;v)=\int \int \int J(\epsilon ,\eta ,w)P(x,y,v| \epsilon ,\eta ,w)d\epsilon \,d\eta \,{dw}.\end{eqnarray} \tag{ 2 }$$

Converting integration variables to polar coordinates ($\epsilon =R\cos \theta$ and $\eta \,=\,R\sin \theta \cos \,i$), Equation (2) can be written as

$$\begin{eqnarray}&&\begin{array}{l}I(x,y;v)\\ =\,\int \int \int J(\epsilon ,\eta ,w)P(x,y,v| R\cos \theta ,R\sin \theta \cos i,{v}_{R,\theta })R\,{dR}\,d\theta .\end{array}\end{eqnarray} \tag{ 3 }$$

The instrumental PSF consists of two terms, the spatial PSF P_s and the velocity spread function P_v. The spatial component of the instrumental PSF is written as

$$\begin{eqnarray}&&{P}_{s}(x,y| \epsilon ,\eta )=\displaystyle \frac{1}{2\pi {\sigma }_{s}^{2}}\exp \left[\displaystyle \frac{{\left(x-\epsilon \right)}^{2}+{\left(y-\eta \right)}^{2}}{2{\sigma }_{s}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

where ${\sigma }_{s}^{2}$ is computed taking into account the FWHM of the minor and major axis of the observed PSF, as well as its positions angle in the sky. The velocity-broadening PSF takes the form

$$\begin{eqnarray}&&{P}_{v}(v| w)=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{v}}}\exp \left[\displaystyle \frac{{\left(v-w\right)}^{2}}{2{\sigma }_{v}}\right].\end{eqnarray} \tag{ 5 }$$

The dispersion of the velocity PSF, σ_v, is given by the thermal and non-thermal line widths added in quadrature (e.g., Hughes et al. 2011),

$$\begin{eqnarray}&&{\sigma }_{v}(r)=\sqrt{\displaystyle \frac{2{k}_{{\rm{B}}}T(r)}{m}+{\xi }^{2}},\end{eqnarray} \tag{ 6 }$$

where k_B is Boltzmann's constant, T(r) is the local disk temperature, m is the average mass per particle, and ξ is a velocity-broadening term adjusted to match the instrumental spectral resolution. No turbulence has been assumed. The local disk temperature T(r) is computed assuming that ${T}_{\mathrm{gas}}(r)\,={T}_{\mathrm{dust}}(r)$ and that the dust grains are in thermal equilibrium with the stellar radiation, i.e.,

$$\begin{eqnarray}&&T(r)={\left(\displaystyle \frac{{L}_{* }}{16\pi \sigma {R}^{2}}\right)}^{1/4},\end{eqnarray} \tag{ 7 }$$

where L_* is the stellar luminosity and σ is the Stefan–Boltzmann constant. The velocity broadening due to Keplerian rotation is estimated by assuming a geometrically flat, azimuthally symmetric circumstellar disk viewed at inclination angle i, in polar coordinates

$$\begin{eqnarray}&&{v}_{r,\theta }=\sqrt{\displaystyle \frac{{{GM}}_{* }}{R}}\cos \theta \sin i.\end{eqnarray} \tag{ 8 }$$

The double-convolution kernel ${\rm{\Pi }}({x}_{k},{y}_{k};{v}_{l}| {R}_{j})$ is defined as the convolution between the spatial PSF P_s and the velocity spread function P_v (Scoville et al. 1983),

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Pi }}(x,y;v| R)\\ =\,\tfrac{1}{2\pi }{\int }_{0}^{2\pi }{P}_{s}(x,y| R\cos \theta ,R\sin \theta \cos i){P}_{v}(v| {v}_{R,\theta })d\theta .\end{array}\end{eqnarray} \tag{ 9 }$$

By combining (1), (3), and (9), the predicted intensity profile at any given position can be calculated from

$$\begin{eqnarray}&&{I}^{n}({x}_{k},{y}_{k};{v}_{l})=2\pi \displaystyle \sum _{J=1}^{{n}_{R}}{R}_{J}{\rho }^{n}({R}_{J}){\rm{\Pi }}({x}_{k},{y}_{k};{v}_{l}| {R}_{j}){\rm{\Delta }}R.\end{eqnarray} \tag{ 10 }$$

Note that the ${I}^{n}({x}_{k},{y}_{k};{v}_{l})$ is integrated over an area corresponding to the instrumental beam. Scoville et al. (1983) showed that, for an axisymmetric disk, the de-projected emissivity distribution ρ(r), can be derived on scales much finer than the instrumental spatial resolution by iterating

$$\begin{eqnarray}&&{\rho }^{n+1}(R)={\rho }^{n}(R)\displaystyle \frac{{\displaystyle \sum }_{k=1}^{{n}_{B}}{\displaystyle \sum }_{l=1}^{{n}_{v}}\left[\tfrac{\widetilde{I}({x}_{k},{y}_{k};{v}_{l})}{{I}^{n}({x}_{k},{y}_{k};{v}_{l})}\right]{\rm{\Pi }}({x}_{k},{y}_{k};{v}_{l}| {R}_{j})}{{\displaystyle \sum }_{k=1}^{{n}_{B}}{\displaystyle \sum }_{l=1}^{{n}_{v}}{\rm{\Pi }}({x}_{k},{y}_{k};{v}_{l}| {R}_{j})},\end{eqnarray} \tag{ 11 }$$

where $\widetilde{I}({x}_{k},{y}_{k};{v}_{l})$ is the observed intensity at position (x_k, y_k) and velocity v_l, and ${I}^{n}({x}_{k},{y}_{k};{v}_{l})$ is the theoretical line intensity profile that would be observed given the surface brightness distribution ${\rho }^{n}(r)$. The positions and velocities x_k, y_k, and v_l are sampled at a total of n_B positions and n_v channels, respectively. The radial emissivity distribution ρ(r) is computed at n_R discrete points spaced by ΔR. The separation ΔR can be much finer than the nominal spatial resolution.

Starting with a constant surface brightness distribution ρ⁰(r) = 1, the theoretical line intensity profile is calculated using Equation (10) and ρ⁰ at n_B positions, and the reduced χ² between the observed and theoretical line intensity profiles is computed as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{{\displaystyle \sum }_{j=1}^{{n}_{B}}{\displaystyle \sum }_{v=1}^{{n}_{v}}{\left({\widetilde{I}}_{j,v}-{I}_{j,v}^{n}\right)}^{2}}{{\sigma }^{2}},\end{eqnarray} \tag{ 12 }$$

where σ corresponds to the rms noise per channel in the observations. We adopt a grid that samples every one-third of the beam to prevent over-sampling and to avoid excessive computational time. For all sources the sampled area is limited to a circular region centered at the stellar position of 2 arcsec in radius.

We iterate Equation (11) until the reduced χ² reaches a value less than unity or when the fractional change in the reduced χ² is less than 1% (which was found to be a good criterion for convergence of the parameters). This approach allows us to derive quickly (generally in less than 10 iterations) the surface brightness distribution of the disk that provides the best fit to the models in comparison to other more time-consuming methods (e.g., radiative transfer modeling and/or visibility computation; Isella et al. 2009; Tazzari et al. 2016). In Figure 2 we show an example of the observed $\widetilde{I}(x,y;v)$ and modeled ${I}^{n}(x,y;v)$ line profiles for the best-fit model for HD 138813 (see Section 3.3).

Zoom In Zoom Out Reset image size

In our modeling approach we adopt a Bayesian method to obtain probability distribution functions of the model free parameters. The parameter search is performed using the Python package emcee (Foreman-Mackey et al. 2013), which uses the affine-invariant implementation of the Markov chain Monte Carlo (MCMC) method to run simultaneously several Markov chains to map the posterior probability distribution. By using many walkers, and by proposing new models based on the relative positioning of the walkers, emcee is effective at handling posterior probability density functions (PDFs) with strong degeneracies.

Our disk model is defined by the disk's position angle θ, inclination i, central velocity v_r, as well as the gas surface brightness distribution ρ(r). Two extra parameters, Δα and Δδ, are added to find the exact centroid position of the disk (e.g., Tazzari et al. 2016). Table 2 shows the fixed model parameters for each system.

At each iteration, the position of a given walker in the parameter space (i.e., one set of model parameters) defines a disk model, for which the best-fit surface brightness distribution is found by Lucy–Richardson deconvolution as described in 3.1. The χ² of this model is computed using its surface brightness distribution, and the resultant χ² value is used to compute the PDF. The next iteration will consider the PDF computed by all walkers in order to define the move to the next point in parameter space. In this manner, the walkers interact with each other and the PDF can be sampled quickly and efficiently.

Typically the MCMC chains are left to evolve during the burn-in phase, during which 1000 walkers sample a broad range of the parameter space in order to locate the maximum of the posterior probability. This is achieved in between 100 and 400 steps. After the burn-in phase, 800 to 1000 more iterations are allowed in order to further sample the posterior probability around the global maximum. After the removal of the burn-in steps, the posterior probability provides information on the marginal distributions of the free parameters (i.e., the uncertainties on the derived free parameters). The advantage of this method is that it can be fully parallelized. Using 24 processors in parallel we were able to sample 1200 iterations of the 1000 individual chains in a week of processing time.

The MCMC search through Lucy–Richardson deconvolution model fitting was run for each of the three gas disks. Figure 3 shows staircase plots of the chains obtained for HD 138813 after the MCMC fitting process. The MCMC staircase plots for HD 131835 and HD 156623 are presented in Appendix A. In Table 3 we show the best-fit parameters for each disk. The uncertainties on the individual parameters are estimated by fitting a Gaussian to the marginalized distributions (corresponding to the diagonal panels of the staircase plots). Figure 4 shows the CO surface brightness distribution computed from the best-fit model of each disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows a comparison between the data and model for the HD 138813 disk, together with the resulting residuals in both integrated intensity (moment 0) and velocity dispersion (moment 1) maps. It can be seen from the residual maps that the models provide a reasonable fit to the observations, with no residuals above 5σ level (15% of the peak flux). The data, model, and residual plots for HD 131835 and HD 156623 are presented in Appendix A.

Zoom In Zoom Out Reset image size

HD 131835 is the most studied of the targets within our sample. Moór et al. (2015) was able to model the ¹²CO(3–2) APEX spectra assuming the gas extends from 35 au (obtained from the continuous disk model in Hung et al. 2015b, and consistent with the location of the first inner ring resolved with SPHERE), and varying the outer radius to fit the single-dish data. For their best-fit model, Moór et al. (2015) measure a systemic velocity of 7.2 km s⁻¹ in the local standard of rest (LSR) frame. The systemic velocity we derive is identical to that measured by Moór et al. (2015). This value of the system's velocity is also consistent with the stellar radial velocity measured in the optical by Rebollido et al. (2018). They measure a heliocentric radial velocity of v_Helio = 2.6 ± 1.4 km s⁻¹, which corresponds to 6.4 ± 1.4 km s⁻¹ after converting to v_LSR. This is also consistent with the estimation of the systemic velocity from the [C i] line data by Kral et al. (2018), who derive a heliocentric velocity of 3.57 ± 0.1  km s⁻¹, which corresponds to a v_LSR of 7.37 ± 0.1 km s⁻¹. The radial velocities derived for HD 138813 and HD 156623 are also consistent with radial velocities measured in the optical. Rebollido et al. (2018) find radial velocities of 7.7 ± 2.0 km s⁻¹ and 4.6 ± 1.5 km s⁻¹ for HD 138813 and HD 156623 respectively (v_LSR).

The inclination and position angles adopted in the models from Moór et al. (2015) were fixed to the values derived from Gemini/GPI observations (74° and 50° respectively; Hung et al. 2015b). More recently Feldt et al. (2017) refined the inclination and PA of the near-IR dust disk to 726${\pm }_{-0.6}^{0.5}$ and 603${\pm }_{-0.2}^{0.2}$, respectively. Kral et al. (2018) fit an inclination of ${76.95}_{2.4}^{3.1}$ to the [C i] line data. The inclination angle of 68.6 ± 1.6 we derive is consistent within 1.7σ with the values from the IR images and also from the ALMA C i data (Kral et al. 2018). For comparison, we compute the inclination of the 1.3 mm continuum dust disk by performing a Gaussian fitting of the 1.3mm data. Using the CASA task imfit to fit in the image plane we derive an inclination angle of 69.8 ± 0.7, whereas fitting in the visibility domain we obtain an inclination angle of 66.4 ± 0.2 (using the CASA task uvmodelfit). The scatter seen in the continuum inclination could be indicative of the limitations of the data or that there are systematics affecting the derivation on the inclination angle and thus the formal errors quoted for the inclination of the gas disk could be underestimated. Differences between the inclinations derived by millimeter and scattered light images are not uncommon. Loomis et al. (2017) reported differences of ∼15° in the disk around AA Tau, which could be attributed to a warp in the inner regions of the disk. Discrepancies in the inclinations of the gas and dust disks measured in the millimeter have also been reported, as is the case of β Pictoris in which the CO is more clumpy than the dust and is located 5 au above the midplane more closely aligned with an inner disk and to the orbit of the planet β Pic b (Matrà et al. 2017a). There are no previous measurements of the disk inclinations for the HD 138813 and HD 156623 disks. Higher-resolution images of the HD 131835 system are necessary to investigate the presence of clumps, warps, or other asymmetries in the gas disk.

It is interesting to note that, similar to the inner and outer radii used by Moór et al. (2015) to fit the ¹²CO(3–2) APEX data, our modeling of the ¹²CO(2–1) line derives a surface brightness distribution which is also confined to a region between 50 and *150 au. This is also similar to the extension of the C i disk (40 to 200 au Kral et al. 2018). To quantify how much the assumption of ${T}_{{\mathtt{gas}}}$ = T${}_{{\mathtt{dust}}}$ could affect the determination of the emissivity distribution we experimented using a temperature profile (T(r) ∝ r^−p with p = 0.4) that yields temperatures of ∼100 K in the inner 20 au, similar to those used in the models of Kral et al. (2017). We find that this does not affect the determination of a cavity in the CO, while the derived disk parameters remain consistent within 2σ–3σ.

Figure 6 shows a comparison between the location of the dust and gas disks for all three targets. The inner and outer radii of the dust emission are taken from Lieman-Sifry et al. (2016), while the location of the gas corresponds to the emissivity profiles from Figure 4. The dust and gas ring around HD 131835 are both confined to a ∼100 au ring. In Figure 6 we have also marked the location of the two brightest dust rings seen in the near-IR with SPHERE (named B1 and B2 by Feldt et al. 2017). The two near-IR dust rings roughly coincide with the peak of the CO emission, which could suggest a common origin. Similarly to HD 131835, HD 138813 shows a ring-like structure in dust (Lieman-Sifry et al. 2016), and also in ¹²CO according to our modeling of the ALMA data. The CO surface brightness from HD 156623 is found to be centrally peaked. Since the inner radius of the dust disk was unresolved at 1.3 mm, it is not possible to conclude whether the peak gas emission resides inside a dust cavity or if it is mixed with the dust. Most of the gas (20% to 60% of the peak surface brightness), however, appears to be co-located with the dust.

Zoom In Zoom Out Reset image size

A simple CO mass estimate can be made by assuming that the emission is in LTE and optically thin. The fractional population of the J = 2 level is ${N}_{2}/{N}_{\mathrm{CO}}=(2J+1){e}^{-16.6/T}/Z(T)$ where $Z(T)\sim T/2.762$ is the rotational partition function for CO (Hollenbach & McKee 1979). For a mean line luminosity of ∼10⁻⁹ L_⊙ as observed (a flux of 1 Jy km s⁻¹ from a disk at 100 pc gives a line luminosity ∼2 ×  10⁻⁹ L_⊙), and with a transition probability ${A}_{21}=6.91\times {10}^{-7}$ s⁻¹, we can estimate the total number of CO molecules from the total line luminosity ${L}_{21}={N}_{2}{A}_{21}{\rm{\Delta }}{E}_{21}$ giving a CO mass

$$\begin{eqnarray}&&{M}_{\mathrm{CO}}^{\mathrm{LTE},\tau \lt 1}=1.25\times {10}^{22}\,\left(\displaystyle \frac{{L}_{21}}{{10}^{-9}{L}_{\odot }}\right)\,T\,{e}^{16.6/T}\quad g\end{eqnarray} \tag{ 13 }$$

and a minimum mass ${M}_{\mathrm{CO}}^{\min }$ of 5.6 × 10²³ g or 10⁻⁴ M_⊕ if all the gas is at ∼16.6 K (corresponding to the upper energy level of the CO(2–1) line).

Deviating from any of the assumptions made above on the optical depth, LTE conditions, or temperature would result in a mass higher than M${}_{\mathrm{CO}}^{\min }$ to reproduce the observed emission. From Equation (13), the mass derived is seen to increase for gas temperatures both higher and lower than 16.6K (∼3 M${}_{\mathrm{CO}}^{\min }$ for 5 K and ∼5 M${}_{\mathrm{CO}}^{\min }$ for 200 K).

The fact that we have spatially resolved emission maps allows us to estimate a density from the mass and place constraints on the validity of the LTE assumption. For the two scenarios considered, we assume that H₂ is the main collision partner for primordial gas and that electrons are the main colliders for secondary origin gas (H atoms from photodissociation of water could be abundant, but these collision rates are lower and electrons dominate; e.g., see Matrà et al. 2017a). For thermally populated levels, the collider gas density needs to be higher than the critical density (${n}_{\mathrm{crit}}=A/\gamma \sim {10}^{4}$ cm⁻³ for H₂, Yang et al. 2010, and ∼100 cm⁻³ for e⁻, Dickinson et al. 1977) everywhere in the disk.

From the radial intensity distribution fits (Section 3.3, Figure 4) and from the emission maps, the CO disk radius is ∼200 au. Assuming a constant disk thickness of 10 au (similar to the Kuiper Belt; Jewitt et al. 1996; Trujillo et al. 2001), the disk volume is ∼10⁶ au³.

The mass needed for LTE can be estimated with a simple constant density assumption as ∼n_crit × volume; the primordial origin scenario therefore requires a CO disk mass of ∼2 × 10²³g (for $n(\mathrm{CO})\sim 1.4\times {10}^{-4}n({{\rm{H}}}_{2})$). This limit is similar to that made with the optically thin assumption above, and LTE is likely a valid assumption for the primordial gas disk. The total disk mass in this case (H₂+CO) is expected to be at least ∼0.01 M_⊕. In the case of secondary origin gas, densities are low enough that LTE may be difficult to attain. Emission is therefore likely to be sub-thermal and Equation (13) therefore likely to underestimate the true disk mass (see also Matrà et al. 2015).

As noted in previous work (Kóspál et al. 2013; Moór et al. 2015; Kral et al. 2017), survival of CO against photodissociation could impose stricter constraints on the mass needed to explain the observed line fluxes, especially for the secondary origin scenario. For an ambient interstellar field (e.g., Habing 1968), the lifetime of a CO molecule is ∼120 yr (Visser et al. 2009). A rough order of magnitude estimate of the column densities available for shielding in the two origin scenarios with masses for emission in LTE as described above yield ${N}_{{{\rm{H}}}_{2}}\sim {10}^{18}$ cm⁻² (primordial, ${n}_{{{\rm{H}}}_{2}}\sim {10}^{4}$ cm⁻³) and ${N}_{\mathrm{CO}}\sim {10}^{16}$ cm⁻² (cometary, n_CO ∼ 100 cm⁻³ ); at these column densities UV shielding of CO by H₂ and CO self-shielding are not very efficient (Visser et al. 2009). Primordial gas disk masses need to be ∼10⁴ higher, ∼25 M_⊕, for CO to be shielded and survive for the age of the system. For the secondary origin scenario, the CO is only marginally self-shielded and the estimated CO mass yields lifetimes of ∼2500 yr; but since there is continuous CO production, the mass in this case depends on the replenishment rate with lower disk masses requiring higher replenishment rates.

To summarize, we find that simple estimates using spatial information from the resolved emission maps and typical CO line luminosities yield plausible masses of a few tens of a M_⊕ of H₂ (and ∼few 10⁻³ M_⊕ of CO) for a primordial gas disk and ≳10⁻³ M_⊕ of CO (depending on the CO production rate) for cometary gas to explain the observed emission. While the above estimates are informative, the emission is quite sensitive to the various simplifying assumptions made; the disk temperature structure and CO photochemistry need to be solved for a more accurate determination of the disk mass and to infer implications for the origin scenarios. Moreover, CO is likely optically thick at these estimated densities, affecting the interpreted mass.

We next use thermochemical models (see Appendix B) that solve for gas temperature and chemistry, consider gas line emission with non-LTE radiative transfer, and fit the observed ¹²CO(2–1) line emission to calculate disk masses.

In these models, we assume typical values for all input parameters and only vary the gas mass to match the observed CO emission. We set initial disk elemental abundances as in the interstellar medium (Jenkins 2009). These values are typical of protoplanetary disk gas with the elemental C abundance relative to H equal to 1.4 × 10⁻⁴. The dust mass, grain size and radial extent are from the spectral energy distribution (SED) modeling of Lieman-Sifry et al. (2016). The dust distribution is kept fixed. We adopt double power laws to describe the surface density distribution of gas,

$$\begin{eqnarray}{\rm{\Sigma }}(r)=\left\{\begin{array}{ll}{{\rm{\Sigma }}}_{0}{\left(\tfrac{r}{{R}_{0}}\right)}^{{p}_{1}} & \mathrm{for}\,{R}_{\min }\lt r\lt {R}_{0}\\ {{\rm{\Sigma }}}_{0}{\left(\tfrac{r}{{R}_{0}}\right)}^{{p}_{2}} & \mathrm{for}\,{R}_{0}\lt r\lt {R}_{\max }\end{array}\right.,\end{eqnarray} \tag{ 14 }$$

where Σ₀ is the surface density at R₀. We vary the surface density distribution of gas, and set the local gas/dust mass ratio accordingly. From the Lucy–Richardson modeling, the inclination and radial extent of the gas disk are also constrained and the only free parameters are the radial dependence of the surface density profile and the integrated gas mass of the disk. We note that while Lucy–Richardson modeling gives the emissivity profile, this does not directly correspond to a surface density, especially when the CO emission becomes optically thick. Gas heating, cooling, and chemistry that result from the adopted surface density distribution are all then calculated by the models (for details see Appendix B). The surface density profile (Σ₀, R₀, p₁, p₂) is then varied as described below until the synthetic model line emission profile matches the observed line emission.

We initially fit the integrated flux by varying the total disk mass in the range 0.1–35 M_⊕ in steps of 5 M_⊕, to narrow down the mass range. We then fit the spatially integrated velocity line profile by varying the disk mass in increments of 0.1 M_⊕, and then select the best-fit models for a more detailed analysis as follows. The surface density exponents p₁ and p₂ are varied manually between −5 and +5 using different step sizes, and then fine-tuned using steps of 0.05. The model data are used to generate synthetic CO line emission data cubes using LIME (Brinch & Hogerheijde 2010), with non-LTE radiative transfer and considering o-H₂, p-H₂, H, and e⁻ as collision partners (Dickinson et al. 1977; Yang et al. 2010; Walker et al. 2015, respectively). The synthetic data cubes are then processed through CASA⁷ version 5.1.0 (McMullin et al. 2007) to compare emission in each velocity channel. The model images are processed through the CASA task simobserve to create synthetic visibilities, using the same integration time, spectral setup, and antenna configuration used for the observations as well as injecting the appropriate amount of thermal noise. The resulting model visibilities are then imaged using the same CLEAN parameters used for the real data.

The integrated spectra for the best-fit models are compared to the ¹²CO(2–1) data in Figure 7. We obtain total H₂ disk masses of 8.1, 7.2 and 10.4 M_⊕, and total CO masses of 5.0 × 10⁻³, 4.2 × 10⁻³ and 3.8 × 10⁻³ M_⊕ for HD 131835, HD 138813, and HD 156623 respectively. The simulated channel maps for the primordial origin model for HD 138813 are compared to the data in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Heating, cooling, and photochemistry in the disk around HD 138813 are described in more detail below; the other two disks are similar in their chemical and physical structure and not discussed. Heating mechanisms considered include collisions with dust, X-rays, and cosmic rays, UV grain heating, H₂ vibrational heating, H₂ formation heating, exothermic chemical reactions, and photoreactions such as the ionization of carbon; cooling is by dust collisions, and several ionic, atomic, and molecular lines. Since disk mass is being inferred using only the CO J = 2–1 line, we need only concern ourselves with heating in the regions where CO resides. None of the three stars have any detected X-rays; we therefore assume an X-ray luminosity of 10²⁷ erg s⁻¹, typical of A stars (e.g., Feigelson et al. 2011). At this level, X-rays do not contribute significantly to the heating. Dust collisions are important in the regions with dust (73–161 au for HD 138813). At typical model gas densities, drag may be sufficient to retain dust grains smaller than the blow-out size (see Appendix C), but we do not consider any small dust population that can be retained to heat the gas (the blow-out size for each disk was taken from Lieman-Sifry et al. 2016). We assume that a small fraction (<1%) of the grains are spatially co-located with the gas outside the main dust belt (this may be possible because of gas drag). We note that such small amounts of dust are still consistent with the SED fitting used to determine dust parameters. Photoelectric grain heating (due to far-UV (FUV) photons) in the absence of very small grains and polycyclic aromatic hydrocarbons (PAHs) is not very significant, and here dust heating is only due to relatively large ∼1μm grains (Kamp & van Zadelhoff 2001). The other main heating mechanisms are UV pumping, H₂ formation, and cosmic rays (we use ${\zeta }_{\mathrm{CR}}\sim {10}^{-16}$ s⁻¹). The cosmic ray rate in the Galaxy could be as high as 10⁻¹⁵ s⁻¹ (e.g., Indriolo & McCall 2012; Neufeld & Wolfire 2017). If we assume ${\zeta }_{\mathrm{CR}}\gtrsim {10}^{-16}$ s⁻¹, then cosmic ray heating typically dominates gas heating. In this case, cosmic rays provide the necessary heating (to excite the CO line) and there is no need to assume that there is any small dust that gets dragged along with the gas. Figure 9 shows the dominant heating in the disk multiplied by the CO number density to emphasize the regions where CO is present. The main coolants in the disk are CO rotational emission and [O i] fine-structure emission. The resulting temperature and density structure is shown in Figure 10. The gas surface density parameters of the best-fit models are presented in Table 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Heating by dust in the 67–147 au region results in an increase in the temperature and therefore higher pressure, leading to a more vertically extended disk. Outside the dust belt, the gas temperature is quite low, and in general insufficient to excite the J = 2–1 line of CO in the outer disk.

Since the primordial origin disk is optically thick in CO, the resulting emission is relatively insensitive to mass and depends on the disk temperature structure. Uncertainties in disk heating mechanisms translate to the calculated gas temperature, and could potentially result in different disk mass estimates for the same line emission flux. However, the disk mass also affects the chemistry and the amount of CO in the disk. H₂, CH, C, and CO have significant cross-sections in the 11.3–13.6 eV energy range and absorb incident UV flux to shield and preserve CO deeper in the disk, with H₂ shielding being the most effective. As discussed earlier, the column densities required are such that the implied masses for CO to survive are high. An additional complication is the paucity of dust in the disk. We assumed that there was dust at the 1% level distributed through the disk; not only does this provide some additional heating to raise the gas temperature outside the dust belt, but also enables formation of molecular hydrogen that can shield CO. Increasing this dust fraction affects the SED-fitting, while lowering it increases the required mass due to inefficient H₂ formation (and CO-shielding). As discussed previously, increasing the cosmic ray ionization rate could also lead to more gas heating. Therefore, although we could drive the disk to lower masses with an increase in gas temperature, especially in the outer disk, and still reproduce the observed emission, the limiting mass necessary for self-shielding makes it very unlikely that the primordial disk mass can be significantly lower than calculated from the models.

Disk masses higher than derived from the modeling are not ruled out. Decreasing the heating in the disk, for example by lowering the cosmic ray rate and/or reducing the dust density, lowers the gas temperature (note that increasing the gas mass reduces the effective cross-section of dust per H nucleus and decreases the gas temperature even as the dust mass is unchanged). Lower temperatures in a more massive disk are consistent with the observed emission as well, and this degeneracy cannot be broken with only one CO transition observed.⁸ However, disk masses cannot be considerably higher than derived here because then the densities become high enough for gas drag to prevent radiation pressure from removing grains in the debris collisional cascade, (see Appendix C). Lack of removal of grains (which is size-independent as both forces depend on grain area) will result in an accumulation of small grains, making the dust disk optically thick and inconsistent with the debris disk classification.

In the cometary model, CO may be released by thermal and UV desorption of ices from dust grains (Grigorieva et al. 2007), or by high-velocity impacts of larger bodies (Zuckerman & Song 2012). In the latter case, the heat generated by the impact would release CO into the gas phase, and the initial temperature of the CO gas would be set by the energetics of the collision that heat the gas. If CO is sublimated from dust, the gas that is released is expected to be approximately at the temperature of the dust grains. After desorption, the gas will initially expand and cool adiabatically and radiatively and will be heated further out in the flow by photo and chemical processes (e.g., Rodgers & Charnley 2002). A full and accurate treatment of the gas temperature requires a chemo-dynamical model, possibly including impact collision modeling, and is beyond the scope of this work. Here, for simplicity, we assume that the gas temperature is equal to the local dust temperature. However, we fully solve for the chemical evolution of the gas; we assume that at the start of the simulation CO and H₂O are released in a 1:10 ratio (Mumma & Charnley 2011), and our initial abundances are such that all the gas is in CO and H₂O. We consider different initial surface density distributions of the gas and solve for time-dependent chemical evolution of the disk after the CO and H₂O are released. We use the same chemical network as described in Appendix B, but now only solve for the dust temperature and set the gas temperature equal to it. For low initial surface densities, molecules are rapidly photodissociated unless there is a sufficient self-shielding. The disk surface density at t = 0 is increased to a point where a column of CO is built up that can explain the observed CO emission. We typically run the chemical model up to times of  2 Myr. The best-fit cometary gas masses are 4.6 × 10⁻³, 3.1 × 10⁻³ and 2.45 × 10⁻³ M_⊕ for CO alone (with total gas content being ∼15 higher) for the disks around HD 131835, HD 138813, and HD 156623. The CO masses derived are higher than the minimum mass optically thin estimate, and almost meet the mass needed for LTE. The temperature of the gas is higher than ∼20 K through most of the disk, and in the outer disk, densities are low enough that non-LTE conditions (collisions with neutrals and electrons) result in sub-thermal emission. The gas surface density parameters of the best-fit models are presented in Table 4.

Figure 11 compares the integrated spectra for best-fit models of cometary origin to the ¹²CO(2–1) data (Lieman-Sifry et al. 2016). The models have been used to generate synthetic spectra using LIME, and then through CASA to produce simulated visibilities and channel map images comparable to the ¹²CO(2–1) ALMA data. The simulated channel maps from the secondary origin model for HD 138813 are compared to the data in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As pointed out by previous studies (e.g., Zuckerman & Song 2012; Kóspál et al. 2013; Dent et al. 2014; Kral et al. 2017), simple estimates suggest that CO survival against photodissociation is limited to ∼120 yr; a continuous replenishment of CO by cometary collisions or outgassing is therefore needed to explain the observed emission. For the best-fit disk masses, the CO column densities are high enough (10^16–17 cm⁻² at ∼80 au) for self-shielding of CO and shielding by other C-species to reduce the photodissocation rate by ∼2 orders of magnitude. The CO photodissociation rate in the disk is moreover not constant, and can decrease in time due to an increase in the gas opacity (due to many C-species, and primarily neutral carbon; see also Kral et al. 2018) built up in the disk as CO is destroyed with time (which however also reduces self-shielding), and this can be calculated from the chemical disk model. CO is also re-formed in the gas phase due to the abundance of O-species (from the destruction of H₂O) which further lowers the net CO destruction rate. In order to maintain a column of gas that can shield CO against UV photons, the production rate needs to be equal to the net destruction rate. We estimate this by considering the initial mass of the disk (at t = 0), the timescale t_f taken to reach the final mass of CO needed to explain the emission, and therefore the production rate $\sim (M(t=0)-M({t}_{f}))/{t}_{f}$. Models with higher initial disk masses lead to higher t_f as expected, and give approximately the same production rates. We note that disks with total gas masses (including water released by collisions) of ∼4 × 10⁻² M${}_{\oplus }$ are needed to obtain CO masses ∼3 × 10⁻³ as derived here, and imply a replenishment rate of ∼1.6 × 10¹⁴ g s⁻¹ (8 × 10⁻⁷ M_⊕ yr⁻¹) to retain the minimum mass of CO needed. Given that to capture the true evolution of the secondary gas disk, one needs to consider chemistry with hydrodynamical flow, we do not attempt any more detailed estimates of the production rate.