Molecules with ALMA at Planet-forming Scales (MAPS). IX. Distribution and Properties of the Large Organic Molecules HC₃N, CH₃CN, and c-C₃H₂

The data presented here were collected as part of the ALMA Large Program MAPS ²⁵ (project ID 2018.1.01055.L, co-PIs, K. I. Öberg, Y. Aikawa, E. A. Bergin, V. V. Guzmán, and C. Walsh). Full details of the scientific scope and targets of the program are provided in Öberg et al. (2021). Hence, we provide only a brief overview here and limit our description to information pertinent to the scope of this work.

MAPS targeted five protoplanetary disks—IM Lup, GM Aur, AS 209, HD 163296, and MWC 480—using four spectral settings: two in Band 3 and two in Band 6. Together these settings covered 15 transitions of the large organic molecules HC₃N, CH₃CN, and c-C₃H₂ (see Figure 1). Observations using two antenna configurations (compact and extended) were used to recover both large- and small-scale emission from the disks. Standard calibration routines were initially performed by ALMA staff, supplemented by additional self-calibration to improve the signal-to-noise ratio (see Öberg et al. 2021 for full details).

Several of our target lines are predicted to be relatively weak; hence, we initially processed the measurement sets of the line-containing spectral windows using a matched filter as described in Loomis et al. (2018b) (see Figure 2). Several Keplerian filters were attempted with varying radial extents, from the full spatial extent of the ¹²CO emission, down to very compact filters matched to the emitting radius seen in test images of the target lines. We used the matched filtering response σ_f to define two classes of line detection: tentative (where 3 < σ_f < 5) and robust (where σ_f > 5). This ensured our criteria for detection were not influenced by our choice of imaging parameters. We further discuss the results of the matched filtering in Section 3.1.

Imaging was performed using the tclean task available in CASA version 6 (McMullin et al. 2007). During cleaning, Keplerian masks were created based on the known geometry of the disks to select regions in position–position–velocity space with emission. The residual scaling technique of Jorsater & van Moorsel (1995) was applied to the final image cubes to mitigate artifacts introduced by the combination of data from multiple antenna configurations. Detailed explanations of these steps are provided in Czekala et al. (2021).

In this work, we make use of the fiducial MAPS imaging products. Given the aforementioned weak nature of many of our targeted lines, coupled with the fact that we compare observations across ALMA Bands 3 and 6, all data presented in this paper utilize images with circularized beams of size ∼03, obtained by applying a uv taper during imaging (see Section 6.2 of Czekala et al. 2021). This was found to provide the optimal combination of angular resolution and sensitivity for our analysis. All image cubes were post-processed in a homogeneous manner through our dedicated pipeline to produce integrated intensity maps (zeroth moment; see Figure 1 and Appendix A) and azimuthally averaged radial emission profiles (see Section 3.3 and Law et al. 2021a, for further details).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to robustly recover emission in the image plane from the weakest transitions, we utilize the line shifting and stacking technique available within the GoFish package, making use of the integrated_spectrum() function (Teague 2019). This approach exploits the known geometry and velocity structure of the disk to deproject the rotation profile and combine Doppler shifted emission to a common centroid velocity reference frame. This results in a single disk-integrated spectrum for each transition that we use to measure a disk-integrated line flux (or a corresponding upper limit in the case of nondetections). Velocity ranges for integration are chosen by eye for each transition to include all positive emission in the main core of the line. Uncertainties are calculated on a per channel basis, taking into account decorrelation along the spectral axis (see also Yen et al. 2016). Figure 3 shows examples of the stacking technique for several transitions in MWC 480, and Table 1 reports these measurements for all transitions and targets.

Zoom In Zoom Out Reset image size

We note that the measurement of line flux for the CH₃CN K = 0 and K = 1 transitions is complicated by their separation of only ∼6 km s⁻¹, leading to significant blending at higher velocities. We attempted to mitigate this by use of a sufficiently tight Keplerian mask, but contamination still occurred (as can be seen in Figure 3). Since these lines are expected to produce symmetric profiles, we therefore adopt the approach shown in Figure 3 for these transitions that involves measuring the line flux from the unblended half of the integrated spectrum and assuming symmetry about the rest velocity of the line in order to recover the full line flux.

We exploit the fact that we have multiple transitions of each molecule to empirically extract both disk-integrated and radial-dependent column densities, N_T , and rotational temperatures, T_rot, wherever possible. Our methodology, which is similar to that presented in Loomis et al. (2018a), can be described as follows. Under the assumption of optically thin emission, the surface brightness of line emission I_ν is related to the column density of molecules in the upper level of each transition ${N}_{u}^{\mathrm{thin}}$ via

$$\begin{eqnarray}&&{I}_{\nu }=\displaystyle \frac{{A}_{\mathrm{ul}}{N}_{u}^{\mathrm{thin}}{hc}}{4\pi {\rm{\Delta }}v},\end{eqnarray} \tag{ 1 }$$

in which A_ul is the Einstein A coefficient and Δv is the intrinsic line width (see, e.g., Goldsmith & Langer 1999). In the case of disk-averaged emission, I_ν = S_ν /Ω, where S_ν is the flux density and Ω is the solid angle subtended by the emission, calculated from the radial extent of each line. Substituting for I_ν and inverting this relation gives

$$\begin{eqnarray}&&{N}_{u}^{\mathrm{thin}}=\displaystyle \frac{4\pi {S}_{\nu }{\rm{\Delta }}v}{{A}_{\mathrm{ul}}{\rm{\Omega }}{hc}},\end{eqnarray} \tag{ 2 }$$

where S_ν Δv is the integrated flux density reported for each of the transitions in Table 1.

This can then be related to the total column density of the emission through the Boltzmann distribution,

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{u}}{{g}_{u}}=\displaystyle \frac{{N}_{T}}{Q({T}_{\mathrm{rot}})}{e}^{-{E}_{u}/{k}_{B}{T}_{\mathrm{rot}}},\end{eqnarray} \tag{ 3 }$$

where g_u and E_u are the degeneracy and energy of the upper level, respectively (see Appendix B), and Q(T_rot) is the partition function at the rotational temperature (linearly interpolated in log-space from calculated values tabulated in the Cologne Database for Molecular Spectroscopy (CDMS); see Müller et al. 2005). Taking the logarithm of Equation (3) yields

$$\begin{eqnarray}&&\mathrm{ln}\displaystyle \frac{{N}_{u}^{\mathrm{thin}}}{{g}_{u}}=\mathrm{ln}{N}_{T}-\mathrm{ln}Q({T}_{\mathrm{rot}})-\displaystyle \frac{{E}_{u}}{{k}_{B}{T}_{\mathrm{rot}}},\end{eqnarray} \tag{ 4 }$$

which can form the basis of a linear least-squares regression that derives the rotational temperature, T_rot, and total column density, N_T , from the best-fitting slope and intercept, respectively.

In the case that τ is not negligible, an optical depth correction factor, C_τ , must be applied such that the true level populations become

$$\begin{eqnarray}&&{N}_{u}={N}_{u}^{\mathrm{thin}}{C}_{\tau },\end{eqnarray} \tag{ 5 }$$

where C_τ = τ/(1 − e^−τ ) and Equation (4) can then be expressed as

$$\begin{eqnarray}&&\mathrm{ln}\displaystyle \frac{{N}_{u}^{\mathrm{thin}}}{{g}_{u}}+\mathrm{ln}{C}_{\tau }=\mathrm{ln}{N}_{T}-\mathrm{ln}Q({T}_{\mathrm{rot}})-\displaystyle \frac{{E}_{u}}{{k}_{B}{T}_{\mathrm{rot}}}.\end{eqnarray} \tag{ 6 }$$

The optical depth of individual transitions can be related back to the upper state level populations via

$$\begin{eqnarray}&&{\tau }_{\mathrm{ul}}=\displaystyle \frac{{A}_{\mathrm{ul}}{c}^{3}}{8\pi {\nu }^{3}{\rm{\Delta }}v}{N}_{u}({e}^{h\nu /{k}_{B}{T}_{\mathrm{rot}}}-1),\end{eqnarray} \tag{ 7 }$$

which implies that C_τ can be written as a function of N_u and substituted into Equation (6) in order to construct a likelihood function ${ \mathcal L }({\rm{d}}{\rm{a}}{\rm{t}}{\rm{a}},{N}_{T},{T}_{{\rm{r}}{\rm{o}}{\rm{t}}})$ to be used for χ² minimization. We use the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al. 2013) to fit the observed data using this likelihood function, generating posterior probability distributions that describe the range of possible values of both N_T and T_rot. Best-fitting values and their uncertainties are chosen from the median and 16th–84th percentile of these posterior distributions, respectively. We assume uniform priors spanning ranges of 10⁸ < N_T < 10¹⁵ cm⁻² and 1 < T_rot < 150 K, respectively, which encompasses the typical values expected for these molecules in disks (see, e.g., Bergner et al. 2018).

The intrinsic line width, Δv, can be described as a combination of thermal and turbulent broadening such that

$$\begin{eqnarray}&&{\rm{\Delta }}v=2\sqrt{\mathrm{ln}2}\sqrt{\frac{2{k}_{B}{T}_{\mathrm{ex}}}{{m}_{X}}+{t}_{0}^{2}\frac{{k}_{B}{T}_{\mathrm{ex}}}{\mu {m}_{{\rm{H}}}}},\end{eqnarray} \tag{ 8 }$$

where m_X and m_H are the masses of the molecule and hydrogen, respectively, μ = 2.37 is the assumed mean molecular weight, and t₀ ∼ 0.01 is the assumed contribution to the line width from turbulence (expected to be small in the line-emitting region of protoplanetary disks; see, e.g., Flaherty et al. 2015). Since Δv is a function of temperature, we iterate the fitting procedure assuming T_ex = T_rot until the resulting best-fitting values for column density and rotational temperature converge, which is usually achieved within ≲5 iterations.

For the determination of radially resolved quantities, we perform the same procedure as for the disk-integrated analysis but where S_ν Δv is obtained from the radial profiles in annular bins of width 0075 (one quarter of a beam) and Ω is the solid angle of each annulus. In these cases, the MCMC fitting procedure is performed independently for each bin, and the results are combined to report values of N_T (r), T_rot(r), and τ(r) for each molecule in each disk.

As explained in Section 2, all line-containing spectral windows were analyzed using a matched filter to confirm the detection of each of our targeted species and lines (see Loomis et al. 2018a, for full details). We used a Keplerian filter that was manually varied in radial extent to achieve the highest filter response (σ_f ) for each molecule. We found that a 200 au radius filter was optimal for HC₃N and c-C₃H₂, whereas a more compact filter with a radius of 100 au provided the strongest filter response for CH₃CN, indicating the more compact nature of the emission from this species. The results of the matched filter are shown in Figure 2.

We strongly detect (σ_f ≳ 10) both lines of HC₃N in four out of five of our sources, but it is not detected in IM Lup. We also strongly detect (σ_f > 10) the 7₀₇–6₁₆ and 6₂₅–5₁₄ transitions of c-C₃H₂ in all disks except IM Lup (where there is only a tentative detection of the 7₀₇–6₁₆ transition). The 6₁₅–5₁₄ transition of c-C₃H₂ is detected with σ_f > 5 in three out of five sources (AS 209, HD163296, and MWC 480).

For CH₃CN, we targeted the K-ladder of the J = 6–5 transition in Band 3 and the J = 12–11 transition in Band 6. In both cases, the K = 0 and K = 1 lines are blended. These transitions of the 12–11 line of CH₃CN are well detected in all sources except IM Lup, for which there is a tentative detection of the 12₁–11₁ transition only (σ_f ∼ 3). Further, we detect emission also from the 12₂–11₂ transition in GM Aur, HD 163296, and MWC 480. The emission from CH₃CN in Band 3 is significantly weaker. We detect the 6₀–5₀ and 6₁–5₁ lines only in AS 209, HD 163296, and MWC 480. There is a tentative detection (σ_f ∼ 3) of these transitions as well as the 6₂–5₂ transition in GM Aur.

In summary, we have detected at least two transitions of each of our targeted species in all disks except IM Lup. We confirm the previously reported detections of CH₃CN and HC₃N in HD 163296 and MWC 480 by Öberg et al. (2015) and Bergner et al. (2018), and report new detections of these species in GM Aur and AS 209, as well as a tentative detection of CH₃CN in IM Lup. Further, we report new detections of c-C₃H₂ in three sources (GM Aur, AS 209, and MWC 480) and confirm the previous detection in HD 163296 reported by Qi et al. (2013).

Figure 1 presents integrated intensity (zeroth moment) maps for the brightest transitions of HC₃N, CH₃CN, and c-C₃H₂ in our target disks (the full gallery of all transitions in all disks is shown in Appendix A).

The c-C₃H₂ emission has a clear ring-like morphology in AS 209, HD 163296, and MWC 480 for both transitions, potentially indicative of an association with the outer dust rings apparent in the continuum images (see Figure 1). The zeroth-moment maps show that emission from this species is weaker in GM Aur and only very tentatively present in IM Lup (confirming the matched filter analysis; see Figure 2). HC₃N also presents a ring-like morphology in both the J = 11–10 and 29–28 transitions toward AS 209 and HD 163296, while only the J = 11–10 transition appears ring-like toward MWC 480 (with the J = 29–28 appearing centrally peaked). GM Aur initially appears to be an outlier with centrally peaked HC₃N emission in both transitions, but we note that a ring-like morphology is also observed toward this source when examining higher-resolution data (015; see Law et al. 2021a). In most cases, the J = 11–10 transition of HC₃N emission is similarly extended to that of c-C₃H₂. However, the higher-energy J = 29–28 transition (E_u = 190 K) appears more compact in all disks in which it is well detected. For CH₃CN, lines in the J = 6–5 ladder are significantly weaker than those in the J = 12–11 ladder. In contrast to the emission morphology for c-C₃H₂ and HC₃N, CH₃CN appears to have a clear ring-like morphology only in AS 209, with a more centrally peaked morphology present in the other sources in which it is well detected (although we note the presence of minor dips of emission in the innermost regions of HD 163296 and MWC 480 on scales comparable to the beam).

Radial profiles are particularly powerful at revealing substructure not immediately apparent in the integrated intensity maps, and allow us to quantify the radial extent of the emission. Figure 4 shows the azimuthally averaged radial profiles of emission (see Law et al. 2021a for details on how these are generated) of all lines detected in four of our sources: AS 209, GM Aur, HD 163296, and MWC 480. IM Lup possessed numerous nondetections, and even those lines that were found were only detected at tentative significance, so we exclude this disk from the subsequent radial analysis.

Zoom In Zoom Out Reset image size

Figure 4 confirms both the relatively compact nature of emission from this suite of large organic molecules as well as the ringed morphology present in many sources. HC₃N exhibits either a ringed morphology (in AS 209 and HD 163296) or a centrally compact morphology (GM Aur and MWC 480). For AS 209, the lower-energy transition (J = 11–10; E_u = 28.8 K) is stronger than the higher-energy transition (J = 29–28; E_u = 190 K), by up to a factor of ∼1.5 at their peak positions. For all other sources, the converse is true, with the higher-energy transition between factors of ∼1.5–3 stronger at the peak of emission. In all disks, the HC₃N emission extends out only to the outer edge of the millimeter dust disk. For AS 209, the emission from the lower-energy transition peaks at ∼50 au, which lies between the third and forth millimeter dust rings (Guzmán et al. 2018b). Similarly, both transitions of HC₃N in HD 163296 peak at ∼40 au, which lies between the first and second millimeter dust rings in the high-resolution Band 6 continuum image (Isella et al. 2018). However, there is a large caveat with regard to this comparison, namely that the synthesized beam of the high-resolution continuum images is significantly smaller than that of the lines (003 versus 03).

The radial profiles for CH₃CN are either centrally peaked (GM Aur and HD 163296) or display a broad ring (AS 209 and MWC 480) where the latter peaks at radii between ∼30–50 au depending on transition. In all disks, the emission extent is well within the millimeter dust. In all cases, the 12₀–11₀ and 12₁–11₁ transitions are the strongest (E_u = 68.9 K and 76.0 K, respectively). These two Band 6 transitions are considerably stronger than the equivalent Band 3 transitions (i.e., 6₀–5₀ and 6₁–5₁; E_u = 18.5 K and 25.7 K, respectively) for all disks by factors of 2–3.

Finally, the radial profiles for c-C₃H₂ present a mostly ring-like morphology in all cases. This emission appears to be related to the location of millimeter dust gaps in AS 209, HD 163296, and MWC 480 (at 60 au, 50 au, and 80 au, respectively), in agreement with the analysis of the higher-resolution (015) imaging products (Law et al. 2021a). However, for AS 209 and HD 163296, the c-C₃H₂ emission peaks outside of that for HC₃N, highlighting that any correlations between dust continuum and line emission can be different for different molecules. The c-C₃H₂ emission is confined to within the outer edge of millimeter emission, except for GM Aur, where there is weak extended emission. In all cases, the 7₀₇–6₁₆ transition is the stronger of the two, although this is mainly reflecting the difference in the Einstein A coefficients for the transitions (see Table 4), as these two transitions have very similar upper energy levels. In summary, there is not a one-to-one relation between line emission from a specific species and dust morphology, nor is there a straightforward correlation between the morphology of emission between different species, except that HC₃N and CH₃CN appear broadly similar across each disk. There are several explanations for the presence of ringed emission in molecular lines in protoplanetary disks (now widely observed; see, e.g., Pegues et al. 2020; Garufi et al. 2021), including a drop in column density in the inner regions due to destructive chemical reactions, an increase in opacity in either dust or lines masking emission from deeper layers in the inner disk regions, or a change in excitation conditions (see, e.g., Facchini et al. 2018; van der Marel et al. 2018; Alarcón et al. 2020). The one constant across the large organic species and disks is that all transitions originate on scales either less than, or comparable to, the extent of the the millimeter dust continuum in each disk.

Here, we present the results of the calculations of both the disk-integrated and radially resolved column densities (N_T ) and rotational temperatures (T_rot) using the methods outlined in Section 2.5. As noted in Table 1, our measurements of the c-C₃H₂ 7₀₇–6₁₆ line are blended with the corresponding para transition (7₁₇–6₀₆) at the same frequency (see Spezzano et al. 2012). For the purposes of the subsequent analysis, we correct the disk-integrated and radially resolved measurements of this line assuming equal contributions from each transition and an ortho-to-para ratio of 3, representative of measurements in other protoplanetary disks (see, e.g., Guzmán et al. 2018a; Terwisscha van Scheltinga et al. 2021; Cleeves et al. 2021).

3.4.1. Disk-integrated Analysis

Figure 5 shows the results of the disk-integrated rotational diagram analysis for the four disks across which multiple transitions of HC₃N, CH₃CN, and c-C₃H₂ are detected. For all analysis in which data from multiple observing epochs and bands are used, we assume an additional 10% error to account for uncertainties in flux calibration. All disk-integrated column densities and rotational temperatures are presented in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Median Values of Column Density (N_T ), Rotational Temperature (T_rot), and Maximum Optical Depth (τ_max) from the Fitting Procedure to the Disk-integrated Rotational Diagrams

	HC3N			CH3CN			c-C3H2
	NT	Trot	τmax	NT	Trot	τmax	NT	Trot	τmax
	(cm−2)	(K)		(cm−2)	(K)		(cm−2)	(K)
IM Lup	<5.5 × 1012	⋯	⋯	<6.5 × 1011	⋯	⋯	<7.3 × 1011	⋯	⋯
GM Aur	${1.9}_{-0.4}^{+0.4}\times {10}^{13}$	${54.4}_{-4.3}^{+5.8}$	0.40	${2.1}_{-0.1}^{+0.2}\times {10}^{12}$	${40.8}_{-2.8}^{+3.2}$	0.04	${1.1}_{-0.1}^{+0.1}\times {10}^{12}$	20–48	0.04
AS 209	${2.9}_{-0.5}^{+0.5}\times {10}^{13}$	${31.1}_{-1.0}^{+1.2}$	1.85	${1.7}_{-0.2}^{+0.2}\times {10}^{12}$	${24.8}_{-1.8}^{+2.3}$	0.11	${7.0}_{-0.2}^{+0.3}\times {10}^{12}$	25–37	0.24
HD 163296	${7.3}_{-1.9}^{+2.5}\times {10}^{13}$	${36.9}_{-2.2}^{+3.1}$	2.62	${2.3}_{-0.2}^{+0.2}\times {10}^{12}$	${35.1}_{-1.9}^{+2.3}$	0.07	${2.1}_{-0.1}^{+0.1}\times {10}^{12}$	24–63	0.03
MWC 480	${7.8}_{-2.7}^{+3.9}\times {10}^{13}$	${37.1}_{-2.9}^{+4.5}$	2.83	${3.5}_{-0.2}^{+0.2}\times {10}^{12}$	${48.6}_{-4.6}^{+5.7}$	0.06	${5.7}_{-0.2}^{+0.3}\times {10}^{12}$	27–69	0.14

Note. Uncertainties correspond to the 16th-84th percentile of the posterior probability distribution of the MCMC fitting procedure.

Download table as:  ASCIITypeset image

We note that, for c-C₃H₂, both targeted transitions are very close in upper energy level. This significantly limits the lever arm over which to calculate a gradient and thus rotational temperature. We therefore limit the range of priors for c-C₃H₂ such that T_mid,0 < T_rot < T_atm,0, where T_mid,0 and T_atm,0 are the midplane and atmospheric temperature at a radial distance of 100 au for each disk, derived from fitting multiple CO isotopologues (see Appendix C, Table 5, and Law et al. 2021b).

Where detected, the disk-integrated column of HC₃N ranges from 1.9 × 10¹³ cm⁻² for GM Aur to 7.8 × 10¹³ cm⁻² for MWC 480. The disk-integrated rotational temperatures are relatively constant across AS 209, HD 163296, and MWC 480 at ∼30–37 K, but somewhat higher in GM Aur at 55 K. It is important to note that the J = 11–10 transition in all disks appears to be close to optically thick with τ = 0.4–2.8, which may suggest there is a reservoir of HC₃N in the disks not probed by our observations.

The disk-integrated column of CH₃CN presents a narrow range with all values lying within 1.7–3.5 × 10¹² cm⁻², though a larger range in rotational temperature is seen for the CH₃CN emission across the disks in the sample with T_rot ∼ 25–49 K. The calculated optical depths indicate that CH₃CN emission is optically thin in all disks (though our radially resolved analysis reveals CH₃CN may be optically thick across some limited regions of the disks; see Section 3.4.2).

For c-C₃H₂, there is around a factor of six spread in the disk-integrated column densities ranging from 1.1 × 10¹² cm⁻² in GM Aur to 7.0 × 10¹² cm⁻² in AS 209. We again note that the rotational temperatures are constrained based on the temperature structure of the disks rather than population diagrams. Based on these assumptions, it appears c-C₃H₂ is optically thin across all disks.

The nondetection of HC₃N and only tentative detections of a single CH₃CN and c-C₃H₂ transition toward IM Lup prevents any determination of well-constrained values for rotational temperature or column density for these molecules. However, by performing the disk-integrated analysis where nondetections provide upper limits, we are able to determine the corresponding upper limit to column density for these molecules in the IM Lup disk (see Table 2). The resulting HC₃N column density of ≲5.5 × 10¹² cm⁻² is a factor of 3–15 times lower than values found in the other disks. Similarly, the upper limit for the column of CH₃CN and c-C₃H₂ in IM Lup are factors of 3–5 and 2–8 times lower than for the other disks, respectively.

3.4.2. Radially Resolved Analysis

The disk-integrated analysis demonstrates interesting similarities and differences across the molecules and sources in our small sample. Since the emission is spatially resolved in most cases, this allows us to also conduct a radially resolved excitation analysis to examine radially dependent variations that cannot be probed by a disk integrated analysis. Figure 6 shows the results of this analysis presenting the radial-dependent rotational temperature, column density, and optical depth profiles for all molecules and sources in which emission is well detected.

Zoom In Zoom Out Reset image size

For HC₃N, the radial column densities reveal some structure, with a plateau of higher column densities typically reached in the inner (≲100 au) region (∼10¹⁴ cm⁻²) declining monotonically with radius to ∼10¹³–10¹² cm⁻² beyond ∼50–100 au. The rotational temperature is also relatively constant across all disks at ∼50 K. There is an indication of an increase in rotational temperature within the inner 20–30 au of GM Aur and MWC 480, although derived quantities in these regions may be affected by issues such as beam smearing. In addition, the emission from the J = 11–10 transition of HC₃N appears to be optically thick throughout a significant radial region of each disk, sometimes as far out as 130 au. Hence, it is possible that there is a significant reservoir of HC₃N in these regions of these disks to which our observations are not sensitive, and so our radially resolved column densities are likely lower limits in these regions.

For CH₃CN, there are more differences seen between disks in the radially resolved column densities than for HC₃N. This is despite the disk-integrated analysis suggesting very similar disk-integrated columns across all four sources. In general, the radially resolved column densities are between a factor of 5–10 higher than those derived from the disk-integrated analysis. For GM Aur and MWC 480, the column density monotonically decreases with radius from ∼5 × 10¹³ cm⁻² to a few 10¹² cm⁻² at approximately 100 au. On the other hand, the column density of CH₃CN in AS 209 and HD 163296 has more of a broad ring-like structure. The rotational temperatures of CH₃CN appear to show a similar radial behavior across all sources, remaining relatively constant at 30–40 K within ∼100 au, which is similar to values derived from the disk-integrated analysis. There are hints of a rise to ∼80 K in the outer regions of GM Aur and AS 209, but this may be due to a limited signal-to-noise ratio in these regions. There is also a significant rise in rotational temperature in the inner region of AS 209 to ∼80 K, as expected based on the weak K = 0 and K = 1 lines across this region (although we note that the uncertainties here are large). Our analysis also reveals that the CH₃CN Band 3 lines, in particular the 6₀–5₀ and 6₁–5₁ transitions, approach optical depths of 1 or higher out to 100 au in three of the disks (AS 209, HD 163296, and MWC 480) and in the inner regions (≲20 au) of GM Aur. This indicates that, similarly to HC₃N, our observations may not be sensitive to the bulk of the CH₃CN emitting material.

For the radially resolved analysis of c-C₃H₂, we are again limited in our calculation of a rotational temperature by the small difference in upper energy level of the targeted transitions. We therefore fix the rotational temperature to the one-dimensional power-law profile, $T(r)={T}_{100}\times {\left(\tfrac{r}{100\,\mathrm{au}}\right)}^{-q}$, fitted to the ¹³CO J = 2–1 brightness temperature for each disk (see Table 5 in Appendix C and Law et al. 2021b), that we extrapolate to all radii covered by our c-C₃H₂ emission. In all cases, the emission is optically thin, or moderately optically thick (τ ∼ 0.6 between 70–100 au in AS 209 and MWC 480). The column density for GM Aur is broadly constant at ∼10¹³ cm², while AS 209, HD 163296, and MWC 480 show broad peaks at values of ∼5–8 × 10¹³ cm⁻² and monotonic decreases with radius beyond ∼100 au.

In summary, our radially resolved analysis shows much more variation across the sources than suggested by a simple disk-integrated analysis. Of particular note are the larger peak column densities extracted for CH₃CN and c-C₃H₂, and the high optical depth of the 11–10 and 6–5 transitions of HC₃N and CH₃CN, respectively, indicating that for these molecules the derived column densities in the inner disk may also be lower limits.

We can use the results of our radially resolved rotation diagram analysis to further examine the origin of the emission from the large organic molecules. Law et al. (2021b) have determined two-dimensional temperature structures, T(r, z), for each of our target disks based on fitting brightness temperatures from the mostly optically thick J = 2–1 transitions of ¹²CO, ¹³CO, and C¹⁸O isotopologues, based on the two-layer model of Dullemond et al. (2020). We can use this temperature structure to examine the spatial origin of the emission from these molecules in each of our disks. Under the assumption that the rotational temperature of each molecule is a measure of its true temperature (a sensible approximation given the high densities in protoplanetary disks), we can use the range of derived T_rot for each molecule from Figure 6 to map the range of heights in the disk from which the emission would originate. We can determine the radial range of emission following the approach of Law et al. (2021a) by calculating the radius encompassing 90% of flux in the radial profile from each molecule. Here, we measure this using the 03 data for consistency with our previous analysis, but these only show minor differences with the results of Law et al. (2021a). Combining these radial (r) and height (z) bounds, we can determine the two-dimensional origin of the emission for each molecule, which is shown for each disk in Figure 7.

Zoom In Zoom Out Reset image size

We note that the original T(r, z) structure derived in Law et al. (2021b) is not a sensitive probe of the midplane temperature structure of the disks. This is because the CO emission from which these temperature structures are calculated originates from intermediate-to-high relative heights in the disks (z/r ≳ 0.1), and temperatures below 20 K were excluded from the fitting procedure. We have therefore renormalized the T(r, z) structures such that the CO snowline (assumed to be 20 K) is reached at a midplane radius consistent with Zhang et al. (2021), whose thermochemical models reproduce both the SED and radial CO column density of each disk. We also note that we exclude the c-C₃H₂ emission from this comparison, since the rotational temperature we assume is not independent of the temperature structure derived from CO in each disk.

Our analysis of the origin of the line emission in T(r, z) space reveals common trends across the four disks. In all cases, the HC₃N emission appears to originate from the highest temperature regime of the two molecules (∼30–80 K), and therefore at higher relative heights, generally in the range of z/r = 0.1–0.4. In contrast, the lower measured temperature of CH₃CN means it originates from lower relative heights, generally z/r ≲ 0.1–0.2. The underlying differences in temperature structure between the warmer disks around Herbig stars (HD 163296 and MWC 480) and the cooler disks around T Tauri stars (GM Aur and AS 209) result in differences in the location of the emitting region of the molecules. In the cooler disks, the large organics primarily trace the warm molecular layer, and may only trace the disk midplane interior to 20 au. In the warmer disks, the large organics trace a wider radial range at the midplane interior to 60 au, and the warm molecular layer exterior to this. In particular, CH₃CN appears to originate almost exclusively from z/r < 0.1 in these warm disks. Nevertheless, the general morphology of the region of emission is similar for all disks—a region close to the midplane interior to the CO snowline, and an elevated warm molecular layer beyond it. Figure 4.1 also shows the extent of the 220 GHz dust continuum emission for each disk (see Sierra et al. 2021). This spans the full radial extent of the plots for each disk (with the exception of the inner ∼30 au for GM Aur), demonstrating that the edge of the millimeter dust in each disk is always exterior to edge of the large organic emission. This is in contrast to, e.g., H₂CO, which appears to be radially extended across all disks studied here (see Guzmán et al. 2021).

We can use the radial column density profiles to estimate the amount of each large organic molecule present in the inner 50 au of each disk (i.e., on scales comparable to the solar system) via conversion to a total surface density and integrating over the corresponding area. These values are listed in Table 3 (left). While there does not appear to an obvious trend for the gas mass of each large organic molecule between the disks (e.g., there appears to be no stellar mass dependence), it appears that HC₃N is significantly more abundant than either CH₃CN or c-C₃H₂ in each disk (by factors of 5–15). We can also estimate the total (gas and ice) mass of each of the large organics. Chemical models of these molecules in disks predict gas-to-ice ratios across a range of values, from 10³ in Walsh et al. (2014) to 10⁵ in Ruaud & Gorti (2019), so we conservatively adopt a 1000 to 1 ratio for our conversion. These values are also shown in Table 3 (right), expressed as a percentage of H₂O ice, estimated from predicted abundances in the thermochemical disk models of Zhang et al. (2021). This again demonstrates that HC₃N is most abundant in each disk. AS 209 emerges as a potential outlier, being particularly rich in large organics when compared to the other disks in our small sample by factors of ∼30 for HC₃N and ∼10 for both CH₃CN and c-C₃H₂. However, this may instead be due to a low modeled H₂O mass, and so observational constraints on this abundance would be required to confirm this.

While determining the precise chemical origin of the large organic species will require dedicated models for each disk, we can use our empirically determined quantities to gain insight into general properties. Following Furuya & Aikawa (2014), we can estimate the typical temperature at which we would expect the large organic molecules to sublimate from icy grain surfaces into the gas phase. For gas densities in the range of 10⁶–10¹² cm⁻³, we would expect this sublimation temperature to lie between 80–150 K. ²⁶ It is therefore interesting that our derived rotational temperatures for HC₃N and CH₃CN (where we again note that the assumed rotational temperature of c-C₃H₂ is tied to the CO emission) are well below this value across the full radial extent of each disk. This suggests that, despite their close association with millimeter dust disk, these molecules cannot originate from direct thermal desorption (with the possible exception of the inner ∼20 au of GM Aur and MWC 480 for HC₃N, although uncertainties are large in these regions). Therefore, these species must be released from ices by nonthermal desorption processes (such as those triggered by cosmic rays and/or X-rays), or be formed in the gas phase directly.

We can also compare the column density ratios of simple and complex molecules in each disk to understand the efficiency of conversion from small to large organic species. Figure 8 shows the radially resolved ratio of column densities for HC₃N/HCN, CH₃CN/HCN, and c-C₃H₂/C₂H obtained from a comparison of our results with those of Guzmán et al. (2021). We indicate the regions across which our analysis of HC₃N or CH₃CN emission suggests τ > 1, and thus where the column density ratio may be a lower limit. We also note that the derivation of HCN and C₂H column density is based on fitting the hyperfine transitions, and is not therefore influenced by the (large) optical depth of the main line component for these molecules. For the HC₃N/HCN ratio, all disks exhibit lower values in their inner (≲50 au) regions of approximately 20% or less. Beyond radii at which the optical depth drops below 1 in GM Aur and HD 163296, the measured ratios increase sharply (even given the increase in associated uncertainties). The CH₃CN/HCN ratio broadly follows a similar pattern to that of the HC₃N/HCN, namely lower values (≲5%) in the inner (≲40–80 au) region. The CH₃CN/HCN ratio rises gradually in the outer regions of all disks, reaching values of 20–30% in the outer regions (>50 au) in all disks, where optical depths are below one. It is interesting to note that the general form of both the HC₃N/HCN and CH₃CN/HCN ratios as a function of radius (e.g., the sharp and shallow rise, respectively) are in agreement with the forward modeling retrieval performed by (Bergner et al. 2018, see their Figure 9).

Zoom In Zoom Out Reset image size

In contrast to the nitriles, the c-C₃H₂/C₂H ratio is relatively flat across all radii in each of the disks at ∼5%. Increases of this ratio to levels of 15–20% are seen in the outer (≳100 au) regions of AS 209 and MWC 480. It is highly likely that c-C₃H₂ is a species that is formed purely in the gas phase in the atmospheres of protoplanetary disks similar to what is found in the interstellar medium (e.g., Loison et al. 2017). However, given the denser conditions within protoplanetary disks, a mechanism is required to maintain a source of carbon (chains) in the gas phase to seed such a chemistry, and this could potentially come from the nonthermal desorption of icy hydrocarbon precursors such as CH₄ and C₂H₆. Hence, this molecule is a counterexample to molecules such as CH₃CN for which it is known that both gas-phase and ice-phase chemistry are needed to explain its origin in protoplanetary disks (e.g., Loomis et al. 2018a). The flat profile of c-C₃H₂/C₂H within 100 au suggests that the gas-phase conversion of simple to complex hydrocarbons is relatively insensitive to radial location in these regions, and occurs at a similar rate in each disk in our sample.

Column density ratios of HC₃N/HCN and CH₃CN/HCN are have been calculated from remote observations of comets (while c-C₃H₂ has not been detected), allowing us to compare the relative chemical complexity of these disks with the organic material in the solar system. Biver & Bockelée-Morvan (2019) collate these values from numerous observations and find HC₃N/HCN ranges from ∼1–80%, and CH₃CN/HCN ranges from ∼3–45%, which we also show on Figure 8. The HC₃N/HCN ratio (or lower limit) is consistent with cometary measurements across all of the disks out to ∼100 au, and higher beyond this. The CH₃CN/HCN ratio is mostly consistent with cometary values, but the the inner regions of the AS 209, HD 163296, and MWC 480 disks are somewhat lower than the ratio measured in comets.

The expected formation zone of comets in the solar system is generally thought to be ≲40 au (see Mumma & Charnley 2011). The large-to-small organic ratios on these scales are difficult to probe in our target disks, due to the spatial resolution of our data (∼30–50 au). However, the general picture that emerges from this comparison is that the outer (50–100 au) regions of all disks are consistent with the composition of cometary material. In particular, the warmer HD 163296 and MWC 480 disks would likely have comet formation zones at correspondingly larger radii, and so this can be reconciled with a "scaled-up" picture of the solar system. If comets were to form in the outer regions of these disks (≳50–100 au), then they will attain a similar composition of nitriles to those observed in the solar system.

An obvious outlier in our small sample is the IM Lup disk, which only exhibits tentative detections of CH₃CN and c-C₃H₂ and no detection of HC₃N. This is in contrast to firm detections of the smaller organic molecules discussed previously (Bergner et al. 2019; Guzmán et al. 2021), which demonstrate the precursors of the large organic molecules are at least present in this disk. IM Lup is the youngest star–disk system in our sample at 0.2–1.3 Myr old (Alcalá et al. 2017). If these large organic molecules are primarily formed in situ within the disk, then IM Lup may have not had sufficient time to build up a detectable gas-phase reservoir. This would be in broad agreement with dark cloud chemical models that demonstrate timescales of 10⁵–10⁶ yr are required to reach peak abundances for these species (see, e.g., Agúndez & Wakelam 2013). Alternatively, if the large organic molecules are primarily inherited from the protostellar phase (see, e.g., Drozdovskaya et al. 2016, 2019; Bianchi et al. 2019; Booth et al. 2021), then a short prestellar collapse phase would result in a lower abundance of these species in the disk.

The physical conditions of the IM Lup disk might also explain its lack of complex molecular emission. This disk has been found to be massive, with total disk mass estimates on the order of 0.17 M_⊙ (Cleeves et al. 2016). This density structure results in an optically thick region at Band 6 frequencies inside ∼50 au, which may be responsible for suppression of emission from CO isotopologues. A similar suppression of emission from other molecules should also occur, and if the emission originates from closer to the midplane (compared to CO), then this suppression may extend across a larger radial region. This scenario of flux deficit may explain the relative weakness of the large organic emission in IM Lup when compared to the other disks in our sample, and is in agreement with detailed studies of line emission from massive disks (Evans et al. 2019). Disentangling the relative importance of each of these processes on the resultant line emission from large molecules in IM Lup will require a detailed disk-specific model.

Bergner et al. (2018) observed MWC 480 with ALMA and performed a rotational diagram analysis of large organic molecular emission. For HC₃N, they found a disk-integrated rotational temperature of 49 ± 6 K, similar to our derived value of ${37.1}_{-2.9}^{+4.5}$ K. However, their column density of 5.8 ± 2.8 × 10¹² cm⁻² is ∼13 times lower than the value we calculate here. This difference stems primarily from the differing spatial resolution of our data, with our 03 observations able to place tighter constraints on the radial extent of the HC₃N emission for a disk-integrated analysis. For CH₃CN, Bergner et al. (2018) find a disk-integrated column density of 1.8 ± 0.4 × 10¹² cm⁻², which is similar to our value of ${3.5}_{-0.2}^{+0.2}\times {10}^{12}$ cm⁻², but a rotational temperature of 73 ± 23 K, somewhat higher than our ${48.6}_{-4.6}^{+5.7}$ K. This could be reconciled if our observations are tracing different emitting layers; Bergner et al. (2018) observed transitions with E_u ∼ 150–250 K, while our addition of the Band 3 data allows us access transitions down to E_u = 20 K.

Loomis et al. (2018a) observed several transitions of CH₃CN toward the T Tauri star TW Hya, performing a rotational diagram analysis across transitions spanning E_u ∼ 70–150 K. They derive a disk-integrated column density of 1.45 ±0.2 × 10¹² cm⁻² and rotational temperature of 32 ± 4 K. While the column density is comparable to our values across all disks, the rotational temperature is generally lower (with the exception of AS 209). This may be indicative of CH₃CN originating from a cooler region of the TW Hya disk compared to the objects studied here, or the result of a cooler disk more generally.

Qi et al. (2013) detected c-C₃H₂ toward HD 163296 using ALMA Science Verification data, finding a single-ring structure from ∼30–165 au. They derive a column density of 2.2 ± 0.2 × 10¹² cm⁻² at 100 au, which is very similar to our disk-integrated value. Cleeves et al. (2021) also recently reported a multi-line analysis of c-C₃H₂ toward TW Hya, and through a forward modeling approach found a best-fit column density of 1–3 × 10¹² cm⁻² with disk-integrated rotational temperatures for the ortho and para form of 55 ± 13 K and 43 ± 14 K, respectively. Such temperatures are comparable to our assumed values from ¹³CO J = 2–1 between 50 and 100 au in three of our disks (with GM Aur being ∼10 K cooler), and the column density range is comparable to our disk-integrated values. However, our radially resolved column density values are a factor of ∼10–50 higher.

Our observationally derived quantities can be compared to chemical models of protoplanetary disks. While disk-specific models encompassing the full chemistry required to explain the abundances of these molecules are not yet available, more general disk chemical models are still informative. Walsh et al. (2014) studied the composition of a representative disk around a T Tauri star with a large gas-grain complex chemical network that included both HC₃N and CH₃CN. Their spatial distribution of gaseous HC₃N and CH₃CN is characterized by an inner, warm component reaching down to the midplane along with a population of molecules found at higher relative heights (z/r ≳ 0.2) for larger radii. While the scales differ, such a morphology is in broad agreement with our derived origin of the line emission in Figure 7, though dedicated radiative transfer modeling would be needed to determine those regions of the disk that contribute the most to the emergent flux for the transitions studied here.

We can also compare our derived molecular abundances to their model via our radially resolved N_T profiles. Across our disk sample, peak N_T (HC₃N) are between 10¹³–10¹⁴ cm⁻² within ∼100 au, but are likely lower limits due to the high optical depth discussed above. Such values are up to 500 times higher than the HC₃N column densities seen in the Walsh et al. (2014) model at similar radii, and even ∼100 times higher than those at 10 au (see their Figure 8 and Table 2). A less extreme picture emerges when comparing the N_T (CH₃CN), where our derived values are ∼5–10 times higher than those in the model between 10 and 100 au. Exploring the chemical structure as a function of height at a radius of 300 au, Walsh et al. (2014) also found that models including only gas-phase chemistry result in the largest fractional abundance of HC₃N achieved in the disk atmosphere, whereas models in which grain-surface (ice) chemistry was included achieved the largest fractional abundance of gas-phase CH₃CN. This demonstrates the importance of a complete treatment of both gas- and ice-phase reactions in such models.

More recently, Le Gal et al. (2019) investigated the effect of changes in C/O ratio in a protoplanetary disk on the resulting abundances of nitrile species including HC₃N and CH₃CN. Their peak column densities for both species (that require C/O = 1.0 and the inclusion of grain-surface formation routes) are approximately 10¹² cm⁻², still a factor of 10–100 smaller than the column densities we derive across our disk sample. Wakelam et al. (2019) recently explored the impact on vertically integrated column densities of molecules in protoplanetary disks excluding and including the effects of dust growth and settling on the disk structure and chemistry and also including ice chemistry. Of particular interest for our work is the sensitivity of the column density values and distribution of HC₃N. In their fiducial models (no grain growth nor settling), they achieve column densities ≳10¹² cm⁻² only within the inner 20 au, with a typical column density of ∼10¹¹–10¹² cm⁻² obtained in the outer disk, and increasing with radius. This is similar to the results of Walsh et al. (2014). However, in models where grain growth and settling are included, they obtain an increase on the order of 1 to 2 orders of magnitude in the column density of HC₃N in the outer disk. The model that obtains the highest peak column density is also a model that begins the chemistry with atomic initial abundances, rather than molecular, indicating also the importance of having free carbon available in the gas phase for the formation of carbon-rich molecules with gas-phase formation pathways such as HC₃N, despite the bulk gas possessing elemental ratio of C/O < 1. It is not specifically explained why dust growth and settling boosts the formation of HC₃N; however, it may to be related to the decrease in altitude of the dust photosphere allowing greater penetration of bond-breaking radiation that then drives a rich gas-phase chemistry.

The broad picture that emerges from the above comparisons is that current static chemical models generally underpredict the abundance of HC₃N and CH₃CN by a significant margin. Only in the case of a gas-grain chemical model that also includes dust growth and settling does the column density, specifically of HC₃N, approach the values derived here (Wakelam et al. 2019). It is also noteworthy that the static models presented in Loomis et al. (2018a) are able to reproduce the radial column density of CH₃CN in TW Hya well, but only when both gas-phase and ice-phase chemistry are included, along with a high photodesorption rate.

Beyond an incomplete picture of chemical formation pathways, it is possible that the above disk chemical models are not complete in their description of physical processes that could cause a higher abundance of organic molecules. For example, the radial drift of large (≳ millimeter-sized) grains from the outer regions of the disk can alter the chemical composition of the inner regions, which may manifest in one of two different ways. If a molecule is formed primarily via grain-surface reactions, then the radial drift of grains will lead to an enhancement of gas-phase abundance within the relevant ice line (see, e.g., Booth et al. 2017). However, if the organic molecule is produced more readily in the gas phase, then radial drift may still play a role. The redistribution of bulk carriers of major elements via radial drift will alter elemental ratios in the inner disk (for example, C/O; see, e.g., Piso et al. 2015) that can alter the formation efficiency of molecules such as HC₃N and CH₃CN (as demonstrated by Le Gal et al. 2019). Quantifying the detailed effects of these dust transport processes on the distribution of large molecules would require coupled models that include radial drift of dust and gas-grain chemical kinetics. While such models employing small chemical networks are beginning to emerge (see, e.g., Booth & Ilee 2019), they have not yet been expanded to include the reactions required to track the chemical evolution of larger molecules.

In addition to the radial motions of solids, vertical motions of gas and dust could also act to alter the abundances of molecules in a disk. Semenov & Wiebe (2011) modeled the chemical evolution of a disk including turbulent transport of gas and dust in the form of turbulent diffusion. They find that the column densities of species we study here—HC₃N, CH₃CN, and c-C₃H₂—can be enhanced by up to factors of 10–20 in the case of a fast mixing scenario compared with a static disk. This enhancement occurs when grains with icy mantles are transported from the disk midplane to warmer vertical layers, causing heavy radical species in the ices to become more mobile. This leads to the production of more complex molecules on the grain surfaces, which are then released into the gas. Models presented in Öberg et al. (2015) to aid the interpretation of the first detections of HC₃N and CH₃CN in MWC 480 also explored the impact of vertical mixing driven by turbulence on the column densities of both species. While high levels of turbulence are in tension with current measurements in the line-emitting regions of protoplanetary disks (α ≲ 10⁻³; e.g., Flaherty et al. 2020) and observations of well-settled millimeter dust (Villenave et al. 2020), it was shown that even models with low levels of turbulent mixing better-reproduced the column densities derived from observations. It is also possible that mixing is driven via local mechanisms such as meridional flows generated by the accretion of embedded planets (Teague et al. 2019); however, the impact of this is yet to be investigated in chemical models of protoplanetary disks.

Despite the impressive diagnostic power of the MAPS observations we present here, it is clear that we are operating at the limits of the data. For the weakest transitions (e.g., CH₃CN 6–5) we are primarily limited by sensitivity, and so future observational campaigns aiming to characterize larger molecules in disks should be designed with this in mind. We have shown that the combination of observations for transitions spanning a large range of upper state energy are particularly powerful, and so this should also be maximized. Further characterization of the gas substructure for the brighter transitions (e.g., c-C₃H₂ 7₀₇–6₁₆) would benefit from deeper, higher angular resolution studies, which may elucidate any (anti)correlations with the millimeter dust substructures that are so well characterized in these disks. Our findings also demonstrate that the current state of the art for protoplanetary disk chemical models may require extension if they are to explain the latest observations. Whether this extension requires the addition of further physical or chemical processes remains to be seen. Nevertheless, it is clear that high angular resolution and sensitivity observations are essential in order to further our understanding of the complex chemistry in disks.