Molecules with ALMA at Planet-forming Scales (MAPS). X. Studying Deuteration at High Angular Resolution toward Protoplanetary Disks

The MAPS Large Program (project code 2018.1.01055.L) used four spectral setups: two in ALMA Band 3 at a wavelength of ∼3 mm and two in ALMA Band 6 at ∼1.3 mm (Öberg et al. 2021). Two array configurations were used: a short baseline configuration and a long baseline configuration, resulting in baselines between 15 and 3638 m. The details of the data calibration are described in Öberg et al. (2021).

2.1.1. CLEANing and JvM Correction

2.1.2. Adopted MAPS Imaging Products

In this paper, we use the following lines observed by MAPS: the J = 3−2 transitions of HCN, DCN, and N₂D⁺ in Band 6 and the J = 1−0 transitions of HCN and H¹³CN in Band 3. For each emission line, images with several different beam sizes were produced by MAPS. For our Band 3 data, we use the MAPS fiducial images with a circular 03 beam. For the Band 6 data, we use images tapered to the same circular 03 beam (Czekala et al. 2021) for the weak DCN and N₂D⁺ lines to improve the sensitivity. For the strong HCN 3−2 line in Band 6, we use both the images with 015 and 03 circular beams, as listed in Table 1. In particular, when fitting for the column densities, we use the image cubes with 03 beams in order to have uniform spatial resolution for all the emission lines included in the fit. Since MAPS provided individual HCN and H¹³CN image cubes with spectral axes centered onto individual hyperfine components (Öberg et al. 2021), we used the CASA tasks regrid and imageconcat to combine these cubes into a single cube. The basic properties of the MAPS data cubes we used are listed in Tables 1 (for HCN, DCN, and H¹³CN) and 2 (for N₂D⁺).

Table 1. Overview of Properties of HCN, DCN, and H¹³CN Data Cubes Provided by MAPS and Used in This Study

Line	Rms a	JvM b	Spectral Resolution	Channel Width	Beam Size	Usage c
	(mJy beam−1)		(km s−1)	(km s−1)
IM Lup
HCN 1–0	0.7	0.72	0.24	0.5	030 × 030
HCN 3–2	0.3	0.25	0.16	0.2	015 × 015	MOM0, PROF
HCN 3–2	0.6	0.49	0.16	0.2	030 × 030	FLUX, AASPEC, COLDENS
DCN 3–2	1.1	0.79	0.20	0.2	030 × 030
H13CN 1–0	0.7	0.73	0.49	0.5	030 × 030
GM Aur
HCN 1–0	1.8	1.00	0.24	0.5	031 × 031
HCN 3–2	0.8	0.57	0.16	0.2	015 × 015	MOM0, PROF
HCN 3–2	1.0	0.74	0.16	0.2	030 × 030	FLUX, AASPEC, COLDENS
DCN 3–2	0.7	0.63	0.20	0.2	030 × 030
H13CN 1–0	1.5	1.00	0.49	0.5	033 × 033
AS 209
HCN 1–0	0.8	0.83	0.24	0.5	030 × 029
HCN 3–2	0.5	0.28	0.16	0.2	015 × 015	MOM0, PROF
HCN 3–2	0.8	0.51	0.16	0.2	030 × 030	FLUX, AASPEC, COLDENS
DCN 3–2	0.7	0.58	0.20	0.2	030 × 030
H13CN 1–0	0.8	0.90	0.49	0.5	030 × 030
HD 163296
HCN 1–0	0.8	0.93	0.24	0.5	030 × 030
HCN 3–2	0.4	0.30	0.16	0.2	015 × 015	MOM0, PROF
HCN 3–2	0.7	0.53	0.16	0.2	030 × 030	FLUX, AASPEC, COLDENS
DCN 3–2	0.9	0.72	0.20	0.2	030 × 030
H13CN 1–0	0.7	0.95	0.49	0.5	030 × 030
MWC 480
HCN 1–0	1.7	0.99	0.24	0.5	031 × 031
HCN 3–2	0.8	0.54	0.16	0.2	015 × 015	MOM0, PROF
HCN 3–2	1.0	0.74	0.16	0.2	030 × 030	FLUX, AASPEC, COLDENS
DCN 3–2	0.8	0.73	0.20	0.2	030 × 030
H13CN 1–0	1.5	1.00	0.49	0.5	031 × 031

Notes.

^a Measured in an emission-free region of the non-primary-beam-corrected data cube. ^b Ratio of the CLEAN beam area to the dirty beam area used in the JvM correction. See Appendix A and Czekala et al. (2021) for details. ^c For HCN 3–2, the following flags indicate the usage: FLUX: disk-integrated flux; MOM0: zeroth moment; PROF: radial emission profile; AASPEC: azimuthally averaged spectra; COLDENS: column density calculation. For the moment 0 maps, the radial emission profiles, and the disk-integrated fluxes, we used the non-primary-beam-corrected images, while primary-beam-corrected images are used for all other analyses.

Download table as:  ASCIITypeset image

Table 2. Overview of Properties and Usage of N₂H⁺ and N₂D⁺ Data Cubes

Line	Rms a	JvM b	δvc	Channel Width	Beam Size	Beam PA d	Data Source e	Usage f
	(mJy beam−1)		(km s−1)	(km s−1)		(deg)
IM Lup
N2H+ 3–2	3.6	0.71	0.15	0.15	042 × 035	76	Qi19
N2D+ 3–2	1.6	0.78	0.09	0.2	03 × 03	N/A	Öberg21	FLUX, MOM0, PROF
N2D+ 3–2	1.8	0.78	0.09	0.2	042 × 035	76	Öberg21	AASPEC, COLDENS
GM Aur
N2H+ 3–2	2.4	0.75	0.15	0.15	031 × 020	−2	Qi19	FLUX, MOM0, PROF
N2H+ 3–2	2.4	0.75	0.15	0.15	032 × 032	N/A	Qi19	AASPEC, COLDENS
N2D+ 3–2	1.2	0.73	0.09	0.2	03 × 03	N/A	Öberg21	FLUX, MOM0, PROF
N2D+ 3–2	1.3	0.73	0.09	0.2	032 × 032	N/A	Öberg21	AASPEC, COLDENS
AS 209
N2H+ 3–2	2.1	0.67	0.15	0.15	034 × 024	−74	Qi19	FLUX, MOM0, PROF
N2H+ 3–2	2.1	0.67	0.15	0.15	035 × 035	N/A	Qi19	AASPEC, COLDENS
N2D+ 3–2	1.1	0.58	0.09	0.2	03 × 03	N/A	Öberg21	FLUX, MOM0, PROF
N2D+ 3–2	1.2	0.58	0.09	0.2	035 × 035	N/A	Öberg21	AASPEC, COLDENS
HD 163296
N2H+ 3–2	2.8	0.84	0.26	0.2	045 × 035	−88	Qi15
N2D+ 3–2	1.2	0.69	0.09	0.2	03 × 03	N/A	Öberg21	FLUX, MOM0, PROF
N2D+ 3–2	1.4	0.69	0.09	0.2	045 × 035	−88	Öberg21	AASPEC, COLDENS
MWC 480
N2H+ 3–2	1.8	0.95	1.21	1.05	093 × 048	−29	Loomis20
N2D+ 3–2	1.3	0.74	0.09	0.2	03 × 03	N/A	Öberg21	FLUX, MOM0, PROF
N2D+ 3–2	2.2	0.74	0.09	0.2	093 × 048	−29	Öberg21	AASPEC, COLDENS

Notes.

^a Measured in an emission-free region of the non-primary-beam-corrected data cube. ^b Ratio of the CLEAN beam area to the dirty beam area used in the JvM correction. See Appendix A and Czekala et al. (2021) for details. ^c Spectral resolution. ^d Position angle of the beam major axis, measured north through east. For circular beams, position angle is undefined. ^e References: Qi19: Qi et al. (2019), Öberg21: Öberg et al. (2021), Qi15: Qi et al. (2015), Loomis20: Loomis et al. (2020). ^f For lines with several image cubes, the following flags indicate their usage: FLUX: disk-integrated flux; MOM0: zeroth moment; PROF: radial emission profile; AASPEC: azimuthally averaged spectra; COLDENS: column density calculation. For the moment 0 maps, the radial emission profiles, and the disk-integrated fluxes, we used the non-primary-beam-corrected images, while primary-beam-corrected images are used for all other analyses.

Download table as:  ASCIITypeset image

No transitions of N₂H⁺ were observed by MAPS. Fortunately, the J = 3−2 transition of this molecule has been targeted by previous observations in ALMA Band 7 (λ ∼ 1 mm) for all five of our targets. For IM Lup, GM Aur, and AS 209, we use data presented in Qi et al. (2019, project code 2015.1.00678.S, ∼03–04 resolution). For HD 163296, we use the data by Qi et al. (2015, project code 2012.1.00681.S, ∼05 resolution). Finally, for MWC 480, we use data presented in Loomis et al. (2020, project code 2015.1.00657.S, ∼1'' resolution). The calibration of the data is described in the original papers. Using CASA 5.6, we imaged the visibilities in the same way as for the MAPS data: first, Keplerian masks (i.e., masks that select the regions of the image cube where emission is expected based on the Keplerian rotation of the disk) were constructed following the procedure described in Czekala et al. (2021), taking into account the hyperfine components listed in Table 7. Using the Keplerian masks, the data were imaged with the multiscale CLEAN deconvolver implemented in the CASA task tclean with scales of [0, 5, 15, 25] pixels and a Briggs parameter of 0.5. Finally, to yield the correct flux, the JvM correction was applied to all images following Czekala et al. (2021). For most sources, this process results in N₂H⁺ image cubes with fairly circular beams. The exception is MWC 480, where the beam is 093 × 048. However, the beam sizes of the resulting N₂H⁺ images do not match the beam sizes of the MAPS N₂D⁺ images. Thus, we produced additional N₂H⁺ and N₂D⁺ data cubes with matching beams using the CASA task imsmooth. These additional cubes are used for the derivation of azimuthally averaged spectra and column densities. For IM Lup, HD 163296, and MWC 480, we only had to smooth the N₂D⁺ data. For AS 209 and GM Aur, both the N₂H⁺ and N₂D⁺ data required smoothing to achieve matching beam sizes. Table 2 lists the properties and the usage of the various N₂H⁺ and N₂D⁺ image cubes.

Disk-integrated fluxes of all the emission lines considered in this work are presented in Table 3. These were calculated by integrating the flux within a Keplerian mask. The outer radii of the masks were determined from visual inspection of the radial emission profiles presented in Sections 3.2 and 3.3. If the radial profile was too noisy to determine an outer edge, we looked for guidance from stronger lines: for the H¹³CN 1–0 lines, we used the same outer radii as for HCN 1–0, and for N₂D⁺ 3–2 toward GM Aur, we used the same radius as for N₂H⁺ 3–2. In addition, for N₂D⁺ 3–2 toward HD 163296 and MWC 480, a nonzero inner radius of the mask was set such that the negative emission in the inner region is avoided (see Section 3.3 and Appendix B). HCN 1–0 and 3–2 have a significant fraction of their flux in hyperfine components. For these lines, we constructed the mask as a superpostion of Keplerian masks, each centered at one of the hyperfine components listed in Table 6. Errors were calculated by repeating the flux measurement procedure at off-source positions and taking the standard deviation of the off-source fluxes. A 10% flux calibration error (Remijan et al. 2021) was added in quadrature. If the resulting S/N was smaller than 3, Table 3 reports the 3σ upper limit. To calculate these disk-integrated fluxes and errors, we used the image cubes that are not corrected for the primary beam (see Remijan et al. 2021, Section 3.5, for a discussion of the primary beam correction). This is because for primary beam–corrected images, the noise increases toward the edges of the image, making our approach to estimate the error from off-source positions invalid. However, we verified that the difference to fluxes extracted from primary-beam-corrected images is negligible. Non-primary-beam-corrected images were also used for the moment 0 maps and the radial emission profiles. For all other analyses presented in this paper, we used the primary-beam-corrected images.

The HCN lines as well as DCN are detected in all five targets. DCN is detected for the first time toward GM Aur. H¹³CN 1–0 is not detected in any of the disks, but a matched filter analysis in the uv plane (Loomis et al. 2018) described in Appendix C yields a tentative detection (∼3σ) for GM Aur. Although undetected (or tentatively detected in the case of GM Aur), the H¹³CN 1–0 data are useful to constrain the HCN column density. N₂H⁺ 3–2 is detected in all five sources. N₂D⁺ 3–2 is firmly detected toward IM Lup, AS 209, HD 163296, and MWC 480, with the detections toward IM Lup and MWC 480 reported for the first time. Toward one source, GM Aur, the S/N of the integrated N₂D⁺ 3–2 flux is only 3.4. The matched filter analysis (Appendix C) and the disk-integrated spectrum (Appendix D.2) also do not provide a definitive detection. Thus, we consider N₂D⁺ 3–2 toward GM Aur tentatively detected.

The MAPS collaboration produced zeroth moment maps by applying a Keplerian mask to the data cubes and integrating over the velocity axis. The Keplerian masks take into account the hyperfine structure of the emission lines. These maps were used for all scientific analyses (in particular to derive radial emission profiles). To mitigate arc-like artefacts in these maps and better visualize radial structures, MAPS also produced "hybrid" zeroth moment maps. These were produced by combining a Keplerian mask and a smoothed σ-clip mask (thus their name "hybrid," see Law et al. 2021a). We emphasize that by using a clipping mask, some emission is inevitably lost if the clipping threshold is larger than 0σ. Therefore, the hybrid zeroth moment maps are for presentational purposes only and are not used for any quantitative analysis. Note also that the use of masks implies that the noise level is not constant over the map. This is because in general, each pixel in the zeroth moment map is calculated by integrating over a different number of channels (see Figure 2 in Law et al. 2021a).

In Figures 1 and 2, we show the hybrid zeroth moment maps produced by the MAPS collaboration for all the emission lines. Several values of the σ clip were tested, and the final value was chosen by visual inspection of the maps. The maps for the archival N₂H⁺ 3–2 data were generated in the same way.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the hybrid zeroth moment maps show various substructures, the radial profiles presented in Sections 3.2 and 3.3 show these substructures more clearly. Thus, we concentrate on discussing the radial profiles in the following.

Figure 3 shows deprojected radial emission profiles produced by the MAPS collaboration for the HCN, DCN, and H¹³CN emission lines. These were produced by azimuthally averaging a zeroth moment map that was produced with a Keplerian mask only (i.e., without a σ-clipping mask). We use the profiles derived by averaging over the full azimuth of the zeroth moment map in order to maximize the S/N. The uncertainty in each radial bin is estimated as the standard error on the mean in the annulus over which the emission was averaged (Law et al. 2021a).

Zoom In Zoom Out Reset image size

The HCN 1–0 and 3–2 emission shows varied radial emission morphologies. Centrally peaked emission is seen in MWC 480, while the other disks show a central depression. Various rings and shoulders are also observed. The morphology of HCN is discussed in more detail in Guzmán et al. (2021) and Bergner et al. (2021). For DCN 3–2, we identify three distinct morphologies:

• 1.  
  Two DCN rings for IM Lup, centered at ∼140 au and ∼350 au, respectively. These rings are also seen in the zeroth moment map (Figure 1).
• 2.  
  Centrally peaked DCN emission and a weak outer ring at ∼280 au for GM Aur.
• 3.  
  A single ring, centered at ∼60 au for AS 209, 30 au for HD 163296, and 70 au for MWC 480. The ring in HD 163296 has a shoulder at ∼100 au.

The relation of the outer DCN rings of IM Lup and GM Aur to other deuterated molecules (N₂D⁺ and DCO⁺) is further discussed in Section 5.3.

The DCN structure identified here is generally consistent with the findings by Law et al. (2021a). However, by employing an azimuthal wedge along the disk major axis (instead of the full azimuth) when calculating the radial profile, they find that the shoulder in the HD 163296 profile is in fact a distinct ring at 118 au. Furthermore, by looking at the higher-resolution (015) image, they identify a ring at 16 au for GM Aur, instead of a centrally peaked profile. For a detailed characterization of the substructures (precise ring centers with uncertainties, gap widths, etc.) see Law et al. (2021a).

Taking into account the additional information gained from the high-resolution profiles by Law et al. (2021a) described in the previous paragraph, there is generally good agreement between the structures seen in the HCN and DCN radial profiles. For IM Lup, the inner DCN ring corresponds to a ring in HCN 3–2, while the outer DCN ring corresponds to a shoulder in the HCN 3–2 profile. For GM Aur, both HCN and DCN show a bright inner ring and a faint outer ring. For AS 209, both HCN and DCN are in a ring, while for HD 163296, both species show a double ring. MWC 480 is a notable exception in that HCN and DCN show different morphologies: while HCN is centrally peaked, DCN is in a ring.

There are also a few associations of the HCN and DCN line emission structure with dust substructures (Law et al. 2021a). For example, the inner HCN and DCN rings in IM Lup are associated with a continuum ring-gap structure at ∼125 au. The most prominent feature is the association of the outer DCN rings in IM Lup and GM Aur with the edge of the dust continuum disk, as discussed in Section 5.3. For a detailed analysis of the relationships between the radial structures seen in the MAPS data for HCN, DCN, as well as other molecules and the dust, see Law et al. (2021a).

Figure 4 shows the radial emission profiles of N₂H⁺ 3–2 and N₂D⁺ 3–2. These profiles are also produced by azimuthally averaging a zeroth moment map generated with a Keplerian mask only (i.e., without a σ-clipping mask). The N₂H⁺ 3–2 radial emission profiles of IM Lup, GM Aur, AS 209, and MWC 480 are consistent with the ones presented in Qi et al. (2019) and Loomis et al. (2020).

Zoom In Zoom Out Reset image size

N₂H⁺ and N₂D⁺ show similar ring emission structures, although their detailed structure differs from source to source. Qi et al. (2019) proposed that the morphology of the N₂H⁺ emission reflects the vertical temperature structure of the disk. A narrow ring with extended tenuous emission (GM Aur) is expected for a disk that is vertically isothermal up to substantial heights (z/r ≈ 0.2) above the midplane. Conversely, disks with a vertical temperature gradient above the midplane are expected to present broad rings (IM Lup, AS 209). For HD 163296 and MWC 480, observations with higher S/N and, in the case of MWC 480, higher angular resolution would be useful to firmly distinguish between these two cases.

For the N₂D⁺ emission, the four sources with a clear detection show ring structures. Most interestingly, IM Lup clearly shows a double ring structure peaking at ∼100 and ∼330 au, similar to the rings detected in DCN (Section 3.2) and DCO⁺ (Öberg et al. 2015). N₂H⁺ in IM Lup also shows a subtle shoulder at the location of the outer ring of N₂D⁺. The relation between the structures seen for DCN, N₂D⁺, and DCO⁺ will be discussed more in detail in Section 5.3.

We note the presence of negative N₂D⁺ emission toward the center of the HD 163296 and MWC 480 disks. This is probably due to the difficulty of achieving a precise continuum subtraction for this line, as discussed in Appendix B.

In order to compute radial column density profiles, we will model the azimuthally averaged spectra of equally spaced, deprojected annuli (radial bins). However, at each spatial pixel of a data cube, the spectrum is shifted with respect to the systemic velocity due to the Keplerian rotation of the gas. Therefore, in the azimuthally averaged spectrum, the emission is spread over a broad range of velocities. This results in a suboptimal S/N. Furthermore, the hyperfine structure of HCN and N₂H⁺ is not resolved in such a broad spectrum and cannot be used to constrain the column density. Therefore, we apply the following procedure prior to azimuthally averaging (Teague et al. 2016; Yen et al. 2016; Matrà et al. 2017): we shift each spectrum by the projected Keplerian velocity, which is given by

$$\begin{eqnarray}&&{v}_{\mathrm{Kep}}(r,\theta )=\sqrt{\displaystyle \frac{{{GM}}_{\star }}{r}}\sin i\cos \theta ,\end{eqnarray} \tag{ 4 }$$

where r and θ are the (deprojected) radial and azimuthal coordinates of the disk, M_⋆ is the (dynamical) stellar mass, and i is the inclination (see Öberg et al. 2021, for the adopted values). Here we are assuming that emission originates from the midplane. The center of the disk is assigned to the proper motion-corrected stellar position. This procedure centers each spectrum at the systemic velocity. When averaging azimuthally, the S/N is increased and the hyperfine structure remains spectrally resolved.

The calculation of the error bars of the averaged spectra is described in Appendix D.1. The size of the radial bins was chosen to be equal to half of the mean of the beam major and minor axes.

The west side of the AS 209 disk is known to be affected by foreground cloud contamination (Öberg et al. 2011; Huang et al. 2016; Guzmán et al. 2018), which should be most relevant for Band 3 data. Thus, for the 1–0 transitions of HCN and H¹³CN, we follow Teague et al. (2018a) and extract all spectra from a ±55° wedge that encompasses the uncontaminated eastern side of the disk. The other disks are not affected by cloud contamination and thus are averaged over the full azimuth.

Figures 5 and 6 show examples of extracted spectra. The full gallery of spectra can be found in Appendix D.2. In Figure 5, the hyperfine structure of HCN 1–0 and 3–2 is readily visible, except in the innermost region, where the finite spatial resolution causes considerable broadening of the spectra. Hyperfine structure is also seen for N₂H⁺ in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To derive the radial column density profiles of HCN and DCN, we use the 3–2 transitions of HCN and DCN covered in Band 6 and the 1–0 transitions of HCN and H¹³CN covered in Band 3. Our modeling includes the hyperfine satellite lines, which help to derive reliable column densities even if the main components are optically thick. Table 6 provides an overview of the hyperfine lines used in the analysis. Combining the 3–2 and 1–0 transitions allows us to constrain the gas temperature.

For each radial bin, we fit the azimuthally averaged spectra calculated in Section 3.4 of all four lines simultaneously, using the image cubes with a circular 03 beam. We start by assuming local thermodynamic equilibrium (LTE), that is, the excitation temperature T_ex is the same for all transitions and equals the kinetic gas temperature T_kin. The free parameters are the excitation temperature, the HCN column density, the DCN column density, and additional parameters describing the line width and velocity offsets of the spectra. This results in a total of nine free parameters (see Table 4). We assume that HCN, DCN, and H¹³CN have the same temperature, that is, that they are cospatial. Although DCN might be more concentrated toward the midplane compared with HCN, their theoretically expected spatial distributions (Aikawa et al. 2018) are similar enough to justify this first-order approximation. The general similarity of the radial emission profiles of HCN 3–2 and DCN 3–2 (Figure 3) also supports this assumption (Huang et al. 2017).

Table 4. Free Parameters for the Fitting of Azimuthally Averaged Spectra to Derive HCN and DCN Column Densities

Parameter	Prior Low a	Prior High b	Unit	LTE/Non-LTE c
Tex	10	100	[K]	LTE
${\mathrm{log}}_{10}{N}_{\mathrm{HCN}}$	3	18 (LTE), 16 (non-LTE)	${\mathrm{log}}_{10}([{\mathrm{cm}}^{-2}])$
${\mathrm{log}}_{10}{N}_{\mathrm{DCN}}$	3	18 (LTE), 16 (non-LTE)	${\mathrm{log}}_{10}([{\mathrm{cm}}^{-2}])$
FWHMB3	0.5 d	variable e	[km s−1]
FWHMB6	0.2 d	variable e	[km s−1]
${\rm{\Delta }}{v}_{\mathrm{HCN}\ 1-0}$	−0.3	0.3	[km s−1]
${\rm{\Delta }}{v}_{{\rm{H}}13\mathrm{CN}\ 1-0}$	−0.3	0.3	[km s−1]
${\rm{\Delta }}{v}_{\mathrm{HCN}\ 3-2}$	−0.2	0.2	[km s−1]
${\rm{\Delta }}{v}_{\mathrm{DCN}\ 3-2}$	−0.2	0.2	[km s−1]
Tkin	10	90	[K]	non-LTE
${\mathrm{log}}_{10}{n}_{{{\rm{H}}}_{2}}$	3	10	${\mathrm{log}}_{10}([{\mathrm{cm}}^{-3}])$	non-LTE

Notes.

^a Lower bound of flat prior. ^b Upper bound of flat prior. ^c Marks parameters used exclusively in the LTE or the non-LTE run. ^d Equal to the channel width. ^e Initial value for innermost radial bin is 20 km s⁻¹. Dynamically adjusted as larger and larger radii are fitted (see Figure 29 and text).

Download table as:  ASCIITypeset image

We fix the H¹³CN/HCN ratio to the ISM value of ¹³C/¹²C = 1/68 (Milam et al. 2005). To explore the parameter space, we use the MCMC method implemented in the emcee package (Foreman-Mackey et al. 2013). The details of the model and the fitting procedure are described in Appendix E. A few example fits are shown in Figure 5 for AS 209, and the full gallery of fits is shown in Appendix D.2.

As can be seen in Figure 5, in the innermost region (roughly within one beam FWHM from the disk center, i.e., two radial bins, corresponding to 30–47 au depending on the disk), the lines are strongly broadened by the velocity gradient within the beam and the hyperfine structure is not resolved (see also Figures 18, 19, 20, 21, and 22 in Appendix D.2). We caution that the strong broadening in the inner two radial bins might introduce additional uncertainties for the inferred column densities that are not reflected by our error bars.

Figures 7 and 8 show the derived optical depths, temperatures, and column densities for all five sources. The HCN 3–2 transition is optically thick (τ₀ > 1) in the inner ∼100 au for all sources except IM Lup.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next performed a fit without assuming LTE, that is, the excitation temperatures are not necessarily equal to the kinetic gas temperature. In this case, we fit for the kinetic gas temperature T_kin and the H₂ number density. Details are again given in Appendix E. The results are shown in Figure 9. The H₂ number density is constrained to ≳10⁶ cm⁻³ for all disks for r ≲ 200–400 au. For comparison, the critical density is 2 × 10⁵ cm⁻³ and 4 × 10⁶ cm⁻³ for HCN 1–0 and 3–2, respectively (Dutrey et al. 1997). The fitted H₂ number density is mostly consistent with the midplane H₂ number density from the MAPS reference disk models (Zhang et al. 2021), shown with black dotted lines. The column densities are not significantly different from the LTE case, as can be seen in Figure 8, where the non-LTE results are overplotted for comparison. This suggests that the emission is indeed in LTE. The exception is the region around 500 au of the disks around IM Lup and HD 163296, where non-LTE conditions might prevail (see Section 5.1.1).

Zoom In Zoom Out Reset image size

We derive column densities of N₂H⁺ and N₂D⁺ with the same procedure as in Section 4.1, assuming LTE. However, since we observed only one transition for each molecule, we did not fit the excitation temperature. Theoretical models predict that the N₂H⁺ abundance peaks around the threshold temperature of CO freeze-out (e.g., Aikawa et al. 2015, 2018), which is constrained by observations to ∼20 K (Qi et al. 2011; Schwarz et al. 2016; Pinte et al. 2018). Thus, for simplicity we fixed the excitation temperature of both molecules to 20 K for our fiducial fits. This results in a total of six free parameters that are fitted for each radial bin: the column densities, the velocity offsets, and the parameter describing line broadening for each of the two lines. These parameters and their priors are listed in Table 5. The hyperfine structure is taken into account in the same way as for HCN and DCN. Tables 7 and 8 list the 28 and 25 hyperfine components considered for the fitting of N₂H⁺ and N₂D⁺, respectively. Figure 6 shows example fits for IM Lup for a few selected radial bins. The complete gallery of fits is shown in Figures 23–27. Similar to the HCN and DCN column densities, we caution that the column densities of the two innermost radial bins might be less reliable than those in the outer regions because of strong line broadening (see, e.g., Figures 23 and 24 for extreme examples of broadening of N₂H⁺ 3–2).

Table 5. Free Parameters for the Fitting of Azimuthally Averaged Spectra to Derive N₂H⁺ and N₂D⁺ Column Densities

Parameter	Prior Low a	Prior High b	Unit
${\mathrm{log}}_{10}{N}_{{{\rm{N}}}_{2}{{\rm{H}}}^{+}}$	3	18	${\rm{l}}{\rm{o}}{{\rm{g}}}_{10}([{\rm{c}}{{\rm{m}}}^{-2}])$
${\mathrm{log}}_{10}{N}_{{{\rm{N}}}_{2}{{\rm{D}}}^{+}}$	3	18	${{\rm{l}}{\rm{o}}{\rm{g}}}_{10}([{{\rm{c}}{\rm{m}}}^{-2}])$
${\mathrm{FWHM}}_{{{\rm{N}}}_{2}{{\rm{H}}}^{+}\ 3-2}$	variable c	variable d	[km s−1]
${\mathrm{FWHM}}_{{{\rm{N}}}_{2}{{\rm{D}}}^{+}\ 3-2}$	0.2 c	variable d	[km s−1]
${\rm{\Delta }}{v}_{{{\rm{N}}}_{2}{{\rm{H}}}^{+}\ 3-2}$	−0.2	0.2	[km s−1]
${\rm{\Delta }}{v}_{{{\rm{N}}}_{2}{{\rm{D}}}^{+}\ 3-2}$	−0.2	0.2	[km s−1]

Notes.

^a Lower bound of flat prior. ^b Upper bound of flat prior. ^c Equal to the channel width (see Table 2). ^d Initial value for innermost radial bin is 20 km s⁻¹. Dynamically adjusted as larger and larger radii are fitted (see Figure 29 and text).

Download table as:  ASCIITypeset image

Figure 10 shows the derived optical depth and column density profiles for all five sources. The N₂H⁺ 3–2 transition is optically thick around the emission peak for IM Lup, GM Aur, and AS 209. The column densities of N₂H⁺ are close to 10¹³ cm⁻² for these sources, while they only reach ≲10¹² cm⁻² for the two warmest disks in the sample around the Herbig stars HD 163296 and MWC 480. For all sources, the N₂D⁺ 3–2 transition is optically thin throughout the whole disk. The N₂D⁺ peak column densities exceed 10¹¹ cm⁻² except for GM Aur.

Zoom In Zoom Out Reset image size

In order to check the dependence of the column density on the assumed excitation temperature, we also run fits assuming that the excitation temperature equals the midplane gas temperature T_mid shown in Figure 30 of Appendix E.1. These temperature profiles are extracted from the MAPS reference models (Zhang et al. 2021). As can be seen in Figure 10, there are some areas where the two assumptions about the excitation temperature result in different column densities. The strongest difference is seen for IM Lup, where the N₂H⁺ column density increases by more than an order of magnitude if assuming T_ex = T_mid. This is because for IM Lup, the model midplane temperature is as low as ∼7 K. However, even with such a high column density, the fit actually underpredicts the flux of the main N₂H⁺ component by a factor of ∼2, while at the same time overpredicting the flux of the hyperfine components (see Figure 23). This strongly suggests that the N₂H⁺ excitation temperature is actually ≳10 K. Beyond ∼300 au where the model midplane temperature rises to ∼15 K, the column densities derived for the two excitation temperature choices agree again. Compared to N₂H⁺, the N₂D⁺ optical depth and column density of IM Lup show a less strong increase if we assume T_ex = T_mid.

Similarly, when assuming T_ex = T_mid, the N₂H⁺ column density toward AS 209 is higher by a factor ∼4 around ∼75 au (where T_mid ≈ 11 K; see Figure 30), but again the N₂H⁺ data are not fit well by this model (see Figure 25), in contrast to the fit where T_ex = 20 K. Finally, Figure 10 shows that, when assuming T_ex = T_mid, the N₂H⁺ column density toward GM Aur is increased by up to a factor of ∼3 in the region between ∼50 au and ∼250 au. In this region, T_mid ≈ 12 K. However, in contrast to IM Lup and AS 209, the model assuming T_ex = T_mid fits the N₂H⁺ data as good as the model assuming T_ex = 20 K (see Figure 24). Thus, if one believes that the temperature in the emitting region is indeed as low as 12 K, the N₂H⁺ column density toward GM Aur would be a factor of ∼3 larger. From a theory point of view, this low temperature is not unreasonable; the models by Aikawa et al. (2015; see their Figure 10) show that, while the N₂H⁺ abundance should peak around the CO freeze-out temperature (∼20 K), it can remain substantial even below the N₂ freeze-out temperature of ∼17 K (see also Aikawa et al. 2021).

4.3.1. DCN/HCN

Using the column density profiles of HCN and DCN calculated in Section 4.1, we now derive the radial profile of the deuteration fraction of HCN. Figure 11 (upper panel) shows the radial profiles of the DCN/HCN ratio. The results from the LTE and non-LTE fits agree well. The inferred DCN/HCN ratios range from ∼10⁻³ to ∼10⁻¹. For IM Lup, the outer DCN ring at ∼350 au is more strongly deuterated than the inner ring at ∼100 au. The data also suggest that the deuteration further decreases inward of the inner ring, although not at high significance. For the disk of GM Aur, the outer ring at ∼300 au is more strongly deuterated than the inner ∼100 au of the disk. AS 209 shows a weak decreasing trend of the deuteration toward the inner parts of the disk. Finally, for the two disks around the Herbig stars HD 163296 and MWC 480, the deuteration in the innermost region is reduced by almost two orders of magnitude compared to the deuteration at ∼150 au. Thus, while there is some diversity in the deuteration profiles toward our five targets, generally the outer disk regions are more strongly deuterated in HCN compared with the inner disk. In the inner 100 au of the disk around GM Aur, the deuteration seems to be increasing toward the star instead, although higher S/N and angular resolution data would be needed for confirmation. These trends are further discussed in Section 5.1.4.

Zoom In Zoom Out Reset image size

4.3.2. N₂D⁺/N₂H⁺

5.1.1. Comparing LTE and Non-LTE

5.1.2. Comparison with Previous Observations

5.1.3. Comparison of HCN and DCN Column Densities to Model Predictions

5.1.4. The Radial Variation of the DCN/HCN Ratio

5.1.5. DCN/HCN as a Function of Temperature

In order to study the contributions of the low- and high-temperature pathways to HCN deuteration, we now look at the deuteration fraction versus the temperature estimated from the multiline analysis. In Figure 12, we show the temperature as a function of radius (top) and the DCN/HCN ratio as a function of temperature (bottom). We only include points with a relative error on the temperature smaller than 20%. This is a somewhat arbitrary choice, but we found that a stricter cutoff results in only a few points available for analysis, while a more generous cutoff results in a large number of data points with poorly constrained temperatures.

Zoom In Zoom Out Reset image size

If, for a given radius, the LTE and non-LTE fits give similar results (relative difference smaller than 15% for both the fitted temperature and DCN/HCN ratio), we only consider the LTE fit. Data points within half a beam FWHM from the disk center are excluded, because the additional line broadening makes the inference of column densities and temperatures more uncertain. We further define that the upper or lower bound of the fitted DCN/HCN is unconstrained whenever the 99th percentile is larger than 10 or the 1st percentile is smaller than 10⁻⁴, respectively. We exclude data points where both bounds are unconstrained by this definition and plot upper and lower limits whenever only the lower or upper bound is unconstrained, respectively.

In Section 5.1.1, we noticed that the gas temperatures derived from the LTE and non-LTE fits differ in certain regions of the disks. For the regions beyond 300 au in IM Lup and HD 163296, these differences can be explained by non-LTE conditions, that is, the excitation temperature does not equal the gas temperature. As a consequence, we only consider the gas temperature derived from the non-LTE fits for r > 300 au for IM Lup and HD 163296, but it turns out that no non-LTE data points from those regions pass the temperature error bar cutoff. For the inner regions of AS 209, HD 163296, and MWC 480 where the LTE and non-LTE temperatures disagree as well, non-LTE seems unlikely, making the inferred temperatures uncertain. Thus, we plot data points from those regions (90 to 140 au for AS 209, 0 to 130 au for HD 163296, and 0 to 70 au for MWC 480) with open symbols to indicate the possibility of systematic uncertainty. We also plot in Figure 12 power-law fits to the ¹²CO 2–1 brightness temperature (Law et al. 2021b) for guidance, since the HCN temperature is unlikely to exceed the CO temperature.

The bottom panel of Figure 12 also shows the H₂D⁺/H₃ ⁺ and CH₂D⁺/CH₃ ⁺ ratios calculated by balancing Equations (1) (low-temperature pathway) and (2) (high-temperature pathway), respectively. We assumed a thermalized ortho-to-para ratio of H₂, which is theoretically expected for protoplanetary disks (Aikawa et al. 2018; Furuya et al. 2019). For the CH₂D⁺/CH₃ ⁺ ratio, we show two curves, corresponding to the exothermicities by Roueff et al. (2013) and Nyman & Yu (2019). The H₂D⁺/H₃ ⁺ and CH₂D⁺/CH₃ ⁺ ratios are upper limits on the HCN deuteration by the low- and high- temperature pathways, respectively, since, for example, the destruction of H₂D⁺ by CO or the destruction of CH₂D⁺ by electrons is neglected. ²⁵ We see that for a considerable number of points, the low-temperature pathway alone cannot explain the observed deuteration fraction. Thus, a contribution from the high-temperature pathway is needed.

Naively, we might expect that the deuteration should decrease with increasing temperature. There are suggestions of such a trend with temperature in the bottom panel of Figure 12 for the data points from AS 209 and MWC 480. However, the presence of two deuteration pathways as well as possible systematic uncertainties in the derived temperatures complicates the interpretation. More data points with a well-constrained temperature and DCN/HCN ratio would be necessary to study the temperature dependence of the deuteration in more detail.

5.1.6. Is the Drop of DCN/HCN in the Inner Disks an Artefact Caused by High Dust Optical Depth?

5.1.7. Comparison with the DCN/HCN Ratio in Comets

5.1.8. Comparison with the DCN/HCN Ratio in Earlier Evolutionary Stages

5.2.1. Comparison with Previous Observation of N₂D⁺ and Theoretical Predictions

5.2.2. N₂D⁺ as CO Snow Line Tracer

In this section, we discuss the prospects of using N₂D⁺ as a CO snow line tracer. Determining the CO snow line from observations of CO emission lines is generally challenging, because even outside of the CO snow line, there is a warm molecular layer with CO gas. This can be seen in Figure 13, which shows the CO column density profiles for the MAPS disks determined from C¹⁸O observations (Zhang et al. 2021). The vertical green bars indicate the model estimates for the CO snow line, that is, the radii where the abundances of CO gas and CO ice become the same in the MAPS reference models (Zhang et al. 2021). In the two disks around the Herbig stars, HD 163296 and MWC 480, the CO column density profiles rapidly increase inward across the CO snow lines, tracing CO evaporation (Zhang et al. 2021). However, no clear trend is seen for the other disks, that is, the snow line cannot be determined from the CO column density profile.

Zoom In Zoom Out Reset image size

Instead, the snow line is usually determined from indirect tracers. In particular, N₂H⁺ has been used to trace the CO snow line (Qi et al. 2013, 2015, 2019). Since N₂H⁺ is destroyed by CO (N₂H⁺ + CO → HCO⁺ + N₂), N₂H⁺ is expected to be abundant outside of the CO snow line. Therefore, the inner edge of the N₂H⁺ emission can be used as an estimate of the CO snow line. However, models by van 't Hoff et al. (2017) and Aikawa et al. (2018) predict that N₂H⁺ can, in addition to the midplane region, also form in a surface layer of the disk. This surface layer can make the inference of the CO snow line location considerably more uncertain (van 't Hoff et al. 2017). Thus, N₂D⁺ might be a better snow line tracer, since no surface layer is expected to form (Aikawa et al. 2018). Here we test how well N₂D⁺ traces the CO snow line.

In Figure 13, we plot the column density profiles of N₂H⁺ and N₂D⁺ (remember that the column densities might carry additional systematic uncertainties for r ≲ 50 au due to severe broadening; see Section 4.2) and the CO snow line estimate derived from N₂H⁺ observations (gray-shaded area, Qi et al. 2015, 2019) and from the MAPS reference models (green bar, Zhang et al. 2021). Except for GM Aur, where the S/N of N₂D⁺ is low, the estimates for the CO snow line approximately correspond to the inner edges of the N₂D⁺ column density profiles. Unfortunately, an accurate measurement of the inner edge of N₂D⁺ can be challenging since it is difficult to achieve an accurate continuum subtraction for the N₂D⁺ 3–2 line, which can result in unphysical negative emission in the inner-disk region (see Appendix B). Furthermore, N₂D⁺ emission is generally weaker than N₂H⁺. Still, our results suggest that N₂D⁺ can in principle be used as an alternative or complement to N₂H⁺ for constraining the CO snow line.

In this section, we discuss ring structures seen in the radial emission profiles of the deuterated molecules DCN, N₂D⁺, and DCO⁺. In Figure 14, we show DCN 3–2 and N₂D⁺ 3–2 from this work (03 resolution) together with DCO⁺ 4–3 for GM Aur (PI C. Qi, project code 2015.1.00678.S, 035 resolution) and DCO⁺ 3–2 for the other disks (Huang et al. 2017; Favre et al. 2019, 026–07 resolution). The DCO⁺ emission profiles were generated from the archival data in the same way as for the MAPS data, that is, by azimuthal averaging of a zeroth moment map generated with a Keplerian mask (Law et al. 2021a).

Zoom In Zoom Out Reset image size

The disks around the T Tauri targets IM Lup, GM Aur, and AS 209 all show outer rings (i.e., a ring beyond the most central emission component) in at least two of the three deuterated molecules. IM Lup shows a double-ring structure in all three deuterated molecules. For GM Aur, both DCN and DCO⁺ have an outer ring at ∼300 au. For AS 209, N₂D⁺ and DCO⁺ show an outer ring at ∼150 au. In addition, by employing a deconvolution technique, Flaherty et al. (2017) revealed an outer DCO⁺ ring at 215 au ²⁷ in HD 163296, in good agreement with the shoulder of the N₂D⁺ emission profile (see also Salinas et al. 2018). Interestingly, all of these outer rings coincide with the edge of the millimeter dust continuum, which is indicated by the vertical gray line in Figure 14. There is also a shoulder in the N₂D⁺ profile of MWC 480 just inside the millimeter dust edge.

The MAPS disks are not the only sources where rings of deuterated molecules are observed. A DCN ring at the millimeter dust edge is also seen in the disk around LkCa 15 (Huang et al. 2017). Furthermore, nondeuterated molecules can show similar associations to the dust edge: Law et al. (2021a) find such behavior among some of the MAPS sources for HCN and H₂CO. Carney et al. (2017), Pegues et al. (2020), and Guzmán et al. (2021) also report examples of enhanced H₂CO emission at the dust edge.

Öberg et al. (2015) suggested that the outer DCO⁺ ring in the IM Lup disk results from nonthermal desorption of CO ice by UV photons at the edge of the dust disk where the UV penetration depth increases. The newly supplied CO gas then reacts with H₂D⁺, which is abundant due to the low temperature in these outer regions, to form DCO⁺. In addition, CO might return to the gas phase at these large radii as a result of a thermal inversion: generic disk models by Cleeves (2016) and Facchini et al. (2017) predict that the removal of large grains due to radial drift leads to an increase in dust temperature at the dust disk edge.

Similar to DCO⁺ forming from photodesorbed CO, an N₂D⁺ ring might be formed from photodesorbed N₂ (Öberg et al. 2009) that reacts with H₂D⁺. To get a ring in both N₂D⁺ and DCO⁺, the CO gas phase abundance should be neither too high (otherwise N₂D⁺ gets destroyed) nor too low (otherwise DCO⁺ will not form efficiently, e.g., Pagani et al. 2011). As for DCN, in cold gas it forms from DCNH⁺ + e⁻ $\longrightarrow$ DCN + H. DCNH⁺ is produced by the reaction of HNC with either DCO⁺ or D₃ ⁺ (Willacy 2007). Like the CO required for DCO⁺ formation, the HNC might originate from UV photodesorption from icy grains or alternatively be formed from other species (e.g., HCN) that are photodesorbed. No direct laboratory measurements of HNC or HCN photodesorption are available, but naively we might expect that the photodesorption yield is roughly similar to other molecules such as CO, H₂O, or N₂ for which those measurements exist (Cuppen et al. 2017). A thermal inversion could also prompt the desorption of N₂ and other precursor molecules of the deuterated species discussed here.

Not all deuterated molecules show rings at the dust disk edge in all five disks. For example, there is no DCN ring at the millimeter dust edges of AS 209 and MWC 480. This suggests that the mechanism(s) for forming deuterated molecules at the dust disk edge depend on disk properties.

We recall that the first step for calculating azimuthally averaged spectra is to shift the spectrum at each spatial position of the data cube according to the projected Keplerian velocity (see Section 3.4). For each velocity v_i , the averaged flux is then computed as an azimuthal average over a deprojected ring within channel i of the (shifted) cube. The error _i on the averaged flux is given by

$$\begin{eqnarray}&&{\epsilon }_{i}=\displaystyle \frac{{\sigma }_{i}}{\sqrt{{N}_{\mathrm{indep}}^{i}}},\end{eqnarray} \tag{ D1 }$$

where σ_i is the standard deviation of the pixels that were averaged and ${N}_{\mathrm{indep}}^{i}$ is the number of independent samples represented by those pixels. It is calculated as

$$\begin{eqnarray}&&{N}_{\mathrm{indep}}^{i}=\sum _{k}\displaystyle \frac{1}{{n}_{\mathrm{corr}}^{i,k}},\end{eqnarray} \tag{ D2 }$$

where the index k runs over all pixels included in the average and ${n}_{\mathrm{corr}}^{i,k}$ is the number of pixels correlated with pixel k (where only pixels included in the average are considered). For example, consider the case where all pixels within a beam are correlated. Then ${n}_{\mathrm{corr}}^{i,k}={N}_{\mathrm{beam}}$ where N_beam is the number of pixels per beam. In that case, ${N}_{\mathrm{indep}}^{i}=N/{N}_{\mathrm{beam}}$ with N being the total number of pixels included in the average. In other words, ${N}_{\mathrm{indep}}^{i}$ equals the number of beams included in the average. However, in general ${n}_{\mathrm{corr}}^{i,k}\lt {N}_{\mathrm{beam}}$ because (1) the deprojected ring over which the average runs can be narrower than the beam, so that for a given pixel, some neighboring pixels are not included in the average, even though they are within a beam distance, and (2) the spectral shifting procedure applied to the data cube prior to azimuthally averaging means that pixels within a beam can be uncorrelated if their relative spectral shift exceeds the spectral resolution (see also Section 2.2 in Yen et al. 2016). Therefore, for a given pixel k, we calculate the number of correlated pixels as the number of pixels that satisfy all of the following conditions: they are in the same channel of the shifted data cube, are within the deprojected ring, are within a beam centered onto pixel k, and have a spectral shift relative to pixel k smaller than the spectral resolution.

As an alternative, we also computed the standard deviation of the azimuthally averaged spectrum measured at velocities without line emission, which we shall denote δ. We find that generally, δ ≈ _i . Note that _i depends on velocity, while δ does not. Most of the time, _i is slightly larger than δ. The exception is the innermost disk region, where _i at times is significantly smaller than δ. This points to an important limitation of estimating the error bar with _i : when most averaged pixels are correlated, ${N}_{\mathrm{indep}}^{i}\approx 1$, and thus _i ≈ σ_i . But σ_i is small if most pixels included in the average are correlated, leading to an underestimation of the error bar. This problem only occurs in the innermost disk where all pixels within the deprojected ring under consideration are part of the same beam. To be conservative, we adopt max(_i , δ) as our final error bar.

Figures 18–22 show the azimuthally averaged spectra of HCN 1–0, HCN 3–2, DCN 3–2, and H¹³CN 1–0 as described in Section 3.4, together with LTE and non-LTE model spectra. Figures 23–27 show the corresponding spectra for N₂H⁺ and N₂D⁺, together with model spectra where T_ex = 20 K and T_ex = T_mid. In Figure 28, we show the azimuthally averaged spectra over the whole disk for N₂D⁺ 3–2. The line is clearly detected toward IM Lup, AS 209, HD 163296, and MWC 480 and tentatively detected toward GM Aur. This is consistent with the disk-integrated flux (Section 3.1) and the results of the matched filter analysis (Appendix C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we consider the column density calculation under the LTE assumption where the excitation temperature T_ex is the same for all transitions of a molecule and equals the gas kinetic temperature T_kin. We further assume that the different molecular species are sufficiently cospatial to share the same temperature. Our model has the following free parameters:

• 1.  
  The excitation temperature T_ex.
• 2.  
  The logarithms of the column densities: ${\mathrm{log}}_{10}{N}_{\mathrm{HCN}}$ and ${\mathrm{log}}_{10}{N}_{\mathrm{DCN}}$. The column density of H¹³CN is fixed to the ISM ratio of ¹³C/¹²C = 1/68 of the HCN column density.
• 3.  
  The FWHM of the Gaussian convolution kernel used to model instrumental broadening of the lines (see below): FWHM_B3 (for Band 3 data, i.e., the 1–0 transitions) and FWHM_B6 (for Band 6 data, i.e., the 3–2 transitions).
• 4.  
  For each line, a velocity offset with respect to the systemic velocity: ${\rm{\Delta }}{v}_{\mathrm{HCN}1-0}$, ${\rm{\Delta }}{v}_{\mathrm{HCN}3-2}$, ${\rm{\Delta }}{v}_{\mathrm{DCN}3-2}$, and ${\rm{\Delta }}{v}_{{\rm{H}}13\mathrm{CN}1-0}$.

This sums up to nine free parameters. For given values of the parameters, model spectra (including hyperfine structure) are calculated by following an approach described in Appendix A of Teague & Loomis (2020). These model spectra can then be compared to the data. The intrinsic emission spectrum is given by

$$\begin{eqnarray}&&{I}_{\nu }(v)=({B}_{\nu }({T}_{\mathrm{ex}})-{B}_{\nu }({T}_{\mathrm{CMB}}))\cdot (1-{e}^{-{\tau }_{\nu }}),\end{eqnarray} \tag{ E1 }$$

with B_ν the Planck function and T_CMB = 2.73 K the cosmic microwave background (CMB) temperature (i.e., the background radiation is assumed to be the CMB). The frequency-dependent optical depth τ_ν is given by

$$\begin{eqnarray}&&{\tau }_{\nu }={\tau }_{0}\sum _{i}{r}_{i}\exp \left(-\displaystyle \frac{{\left(v-{v}_{\mathrm{sys}}-\delta {v}_{i}-{\rm{\Delta }}{v}_{X}\right)}^{2}}{2{\sigma }_{{\rm{v}}}^{2}}\right),\end{eqnarray} \tag{ E2 }$$

with v_sys the systemic velocity (taken from Öberg et al. 2021) and Δv_X the free parameter describing the velocity offset of the observed emission line X. The sum goes over all hyperfine components, with r_i the strength and δv_i the velocity offset of each component. The strength r_i is given by

$$\begin{eqnarray}&&{r}_{i}=\displaystyle \frac{{g}_{u,i}{A}_{{ul},i}}{\sum _{j}{g}_{u,j}{A}_{{ul},j}},\end{eqnarray} \tag{ E3 }$$

with g_u the weight of the upper level and A_ul the Einstein coefficient for spontaneous emission of the hyperfine line. Tables 6, 7, and 8 list the values of molecular parameters used in this work. τ₀ is the optical depth of the hyperfine-unresolved transition at the line center, given by

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{{hc}}{4\pi }({x}_{l}{B}_{{lu}}-{x}_{u}{B}_{{ul}})N\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{v}}}},\end{eqnarray} \tag{ E4 }$$

with N the column density. Here the Einstein B coefficients and the fraction of molecules in the upper or lower level (x_u and x_l ) refer to the hyperfine-unresolved transition. The latter are given by

$$\begin{eqnarray}&&{x}_{u}={g}_{u}\displaystyle \frac{{e}^{-{E}_{u}/({{kT}}_{\mathrm{ex}})}}{Z({T}_{\mathrm{ex}})},\end{eqnarray} \tag{ E5 }$$

$$\begin{eqnarray}&&{x}_{l}={g}_{l}\displaystyle \frac{{e}^{-{E}_{l}/({{kT}}_{\mathrm{ex}})}}{Z({T}_{\mathrm{ex}})},\end{eqnarray} \tag{ E6 }$$

where the partition function Z is taken from the CDMS database. ²⁸

In contrast to Teague & Loomis (2020), the width of the intrinsic spectrum is fixed to the thermal width, that is,

$$\begin{eqnarray}&&{\sigma }_{v}=\sqrt{\displaystyle \frac{{{kT}}_{\mathrm{kin}}}{m}},\end{eqnarray} \tag{ E7 }$$

where the kinetic temperature T_kin = T_ex in LTE and m is the mass of the molecule. We ignore turbulent broadening given recent results suggesting low turbulence (e.g., Teague et al. 2018b). Finally, the spectrum given by Equation (E1) is convolved with a Gaussian kernel with FWHM given by the free parameter FWHM_Bx. This models any broadening of the line by effects such as beam smearing or imprecise Keplerian shifting. This broadening is highest in the inner region of the disk and gradually becomes smaller at larger radii (see Figures 18–27).

Table 4 shows an overview of the free parameters. We explore the parameter space using the MCMC method implemented in the emcee package (Foreman-Mackey et al. 2013) for each radial bin separately, starting with the innermost bin. We run 200 walkers, each with a total of 5000 steps. The walkers of radial bin n + 1 are initialized using the result from radial bin n. The first 2500 steps are discarded prior to analyzing the results. We choose flat priors over the ranges detailed in Table 4.

To take into account the correlation between neighboring wavelengths in the azimuthally averaged spectra, we rescale the error bars in the spectra by multiplying by the square root of the correlation length (Booth et al. 2017). We define the latter as max(1, D_corr), where D_corr is the FWHM (in units of spectral pixels) of the autocorrelation of the spectrum. Before calculating the autocorrelation, we remove the region of the spectrum with line emission and then subtract the mean value such that the line-free spectrum has a mean value of zero. The correlation length is typically ≲3 spectral pixels, meaning that the rescaling of the error bars is marginal. However, we find that the correlation length can be much larger (up to ∼20 spectral pixels) for the azimuthally averaged spectra extracted from the innermost beam. This is probably because the individual spectra that were azimuthally averaged are all spatially correlated. This spatial correlation is transformed to spectral correlation by the Keplerian shifting procedure.

The correlation length is difficult to determine for spectra extracted from inner-disk regions because the line emission becomes very broad. Thus, no line-free region is available. To circumvent this problem, for the innermost regions (within 1.5 beams), we adopt the correlation lengths determined from the weakest lines: H¹³CN 1–0 for HCN 1–0, HCN 3–2, and DCN 3–2, and N₂D⁺ 3–2 for N₂H⁺ 3–2.

As we move outward in radius, the upper boundaries of the priors on FWHM_B3 and FWHM_B6 are adjusted according to the following scheme. First, we calculate the S/N of FWHM_Bx, defined as the median of the MCMC ("signal") divided by the mean of the 16th and 84th percentiles ("error"). If the S/N is higher than a set threshold (we choose between 3 and 8), then the upper boundary of the prior for the next radial bin is set to the median plus five times the error. This procedure reflects the observed sharp decrease of the line width with increasing radius and allows us to set tighter priors in the outer regions that show little signal. In the case where only the B6 or B3 FWHM exceeds the S/N required for adjustment, we still adjust both FWHM by using the higher S/N FWHM as a basis to estimate a conservative upper boundary for the lower S/N FWHM, taking into account the difference in channel size if necessary. Figure 29 shows an example of the fitted FWHM together with the radially adjusted prior boundaries.

Zoom In Zoom Out Reset image size

To derive the column densities of N₂H⁺ and N₂D⁺, we used an analogous approach and simultaneously fitted the 3–2 transitions. However, the excitation temperature was fixed to 20 K. In addition, to test the dependence on the assumed temperature, we also run fits with the excitation temperature fixed to the midplane temperature extracted from the MAPS reference models by Zhang et al. (2021; see our Figure 30). These source-specific thermochemical models fit the spectral energy distribution, (sub)millimeter continuum images, and vertical locations of the CO-emitting layers reported by Law et al. (2021b). We note that the increase of temperature in the outer disk (particularly clearly seen for IM Lup at ∼300 au) coincides with the truncation radius of the pebble disk, which makes the midplane warmer.

Since the N₂H⁺ and N₂D⁺ data come from different observations, two separate parameters for the FWHM of the Gaussian kernels were used, one for each line. The free parameters and their priors are listed in Table 5.

The method described here is essentially the same as employed by Bergner et al. (2021) to calculate HCN and CN column densities and by Guzmán et al. (2021) to calculate HCN and C₂H column densities. One difference is that Bergner et al. (2021) and Guzmán et al. (2021) do not convolve the intrinsic spectrum but simply redistribute the total flux of the intrinsic spectrum in a model spectrum with broader components. They also do not include H¹³CN in their fits of the HCN column density. Furthermore, Guzmán et al. (2021) employ a two-step fitting method, where they first fit HCN 1–0 and 3–2 at 03 resolution to constrain the excitation temperature and then fit HCN 3–2 at higher 015 resolution, using the constraints on T_ex from the first step. The HCN column densities derived by Bergner et al. (2021) and Guzmán et al. (2021) are generally consistent with our results, with differences only seen inward of 50–100 au (Guzmán et al. 2021, Appendix C).

We repeated the calculation of the HCN and DCN column densities described in Appendix E.1 without the assumption of LTE, that is, T_ex is not necessarily equal to T_kin and can now take different values for each of the emission lines. We now fit the kinetic gas temperature T_kin and the H₂ number density ${n}_{{{\rm{H}}}_{2}}$ with priors detailed in Table 4. Instead of using Equations (E5) and (E6), the fractional-level populations are now determined by solving the statistical equilibrium of the excitation and deexcitation processes. To this end, we use the radiative transfer code pythonradex. ²⁹ For the purpose of the radiative transfer, we ignore the hyperfine structure: we assume that all hyperfine lines of a given emission line share the same excitation temperature (Teague & Loomis 2020) and we do not consider cross-excitation. The excitation temperature for a given transition is then computed from the fractional-level populations by using

$$\begin{eqnarray}&&\displaystyle \frac{{x}_{u}}{{x}_{l}}=\displaystyle \frac{{g}_{u}}{{g}_{l}}{e}^{-{\rm{\Delta }}E/({{kT}}_{\mathrm{ex}})},\end{eqnarray} \tag{ E8 }$$

where ΔE = E_u − E_l . The rest of the calculation is identical to the LTE case. In practice, we precompute a grid of fractional-level populations and excitation temperatures in the T_kin-${n}_{{{\rm{H}}}_{2}}$-space and interpolate on that grid during the MCMC run for computational efficiency. We adopt the scaled HCN-He collision rates by Dumouchel et al. (2010) available in the LAMDA database ³⁰ for all three species (HCN, DCN, H¹³CN). We find that the calculation of the non-LTE level populations fails for column densities above ∼10¹⁶ cm⁻². Thus, we reduce the upper boundary of the column density priors compared to the LTE case (Table 4). Because we find column densities below 10¹⁶ cm⁻² with both the LTE and non-LTE runs, this should not have a significant effect on our results.