Exploring HNC and HCN line emission as probes of the protoplanetary disk temperature

The main input for physical parameters is the disk density structure. We adopt the gas surface density profile as for a viscous evolving accretion disk (viscosity $\nu\propto R^{\gamma}$, Lynden-Bell & Pringle 1974; Hartmann et al. 1998),

where $\Sigma_{\rm c}$ is the surface density at the characteristic radius $R_{\rm c}$ and is set to yield the required total disk mass. The power-law index $\gamma$ is taken as 1 (Hughes et al., 2008; Andrews et al., 2011) and assumed to not change with time. The gas density in the vertical direction follows a Gaussian distribution, with a scale-height angle of $h=h_{c}(R/R_{c})^{\psi}$. The physical scale height is thus given as $H\sim Rh$. The dust surface density is scaled from the gas surface density assuming a typical gas-to-dust ratio of 100. Because larger grains are more settled toward the midplane, we consider two populations of grains, following D’Alessio et al. (2006): a small-grain population (0.005-1 $\mu$m) with the same scale height $h$ as the gas, and a large-grain population (1-1000 $\mu$m) with a reduced scale height of ${\chi}h$, where ${\chi}=0.2$. The fraction of dust surface density distributed to the large grains is described by $f_{\rm large}$ such that $\Sigma_{\rm dust}=f_{\rm large}\Sigma_{\rm large}+(1-f_{\rm large})\Sigma_{\rm small}$. We have tested models in which we varied the $f_{\rm large}$ parameter between 0.9 and 0.99, but found no significant differences in line emission of cyanides (Cazzoletti et al., 2018). We took $f_{\rm large}=0.99$ in our models because large grains dominate the dust mass budget.

High-energy radiation from the central star and interstellar medium affects the disk chemistry. The stellar spectrum is modeled as a blackbody for a given effective temperature $T_{\rm eff}$ and stellar luminosity $L_{\ast}$. The FUV spectrum (6–13.6 eV) is particularly important for the chemistry of cyanides, therefore we considered two types of stars: a T Tauri star with $T_{\rm eff}$=4000 K and $L_{\ast}$=1 $L_{\odot}$, and a Herbig Ae star with $T_{\rm eff}$=10000 K and $L_{\ast}$=10 $L_{\odot}$, which we refer to as T Tauri and Herbig models throughout the paper. Additional UV flux from accretion can also be included for T Tauri stars, modeled as a blackbody emitting at 10000 K on to a 1$M_{\odot}$ and 1.5$R_{\odot}$ star, with the total accretion luminosity controlled by the mass accretion rate $\dot{M}_{\rm acc}$. We took the typical accretion rate for young T Tauri stars of $10^{-8}\,\rm M_{\odot}\,yr^{-1}$ (Hartmann et al., 2016; Manara et al., 2016) as the default in our model and varied the value to explore the effects of UV radiation on cyanide chemistry (see Figure 1 of Visser et al. 2018 for the UV flux change with varying accretion rates and Figure 11 of Cazzoletti et al. 2018 for the consequences on CN emission). We set an interstellar UV flux of $G_{0}$, as defined in Draine (1978), of $\rm\sim 2.7\times 10^{-3}\,erg\,s^{-1}\,cm^{-2}$ between 6 and 13.6 eV range, and a cosmic-ray ionization rate of $\rm 5\times 10^{-17}\,s^{-1}$ per $\rm H_{2}$. The X-ray spectrum was taken as thermal radiation at $7\times 10^{7}$ K between 1 and 100 keV, with the X-ray luminosity of 10³⁰ erg s^-1. The same prescription was used for the T Tauri and Herbig models. While abundances of some molecules or molecular ions such as HCO⁺ may be sensitive to the choice of X-ray luminosity or temperature, parameter studies have shown that most species, including HCN, HNC, and their ratio, have only small variations over the wide range of X-ray values that is applicable to T Tauri and Herbig stars (Stäuber et al., 2005; Bruderer et al., 2009; Bruderer, 2013; Cleeves et al., 2017). A small grid of models was performed by varying disk mass, disk vertical height, and flaring for the T Tauri and Herbig models to represent a wide range of physical and chemical conditions. The disk and stellar parameters used in the models are listed in Table 1.

Our models used the new reduced network presented in Visser et al. (2018) for CNO chemistry, which is based on the network of Bruderer et al. (2012) and includes the reactions for cyanides discussed in Loison et al. (2014). The entire chemical model is the same as was used for the investigation of CN emission in disks (Cazzoletti et al., 2018), except for two changes made relating to the destruction of HNC: (1) we used an updated reaction rate coefficient of HNC+H based on Graninger et al. (2014) with an energy barrier of 200 K, and (2) we included the reaction of HNC+O with an energy barrier of 20 K based on the recent study of the HNC-to-HCN ratio in molecular clouds by Hacar et al. (2020). The updated network has negligible effects on the CN emission. Reactions involving isotopologs were not considered. The same photodissociation rate was used for HCN and HNC; the actural HNC photodissociation rate may be a factor of 2 higher than that of HCN (Aguado et al., 2017).

The network contains gas-phase and grain-surface reactions, including neutral-neutral and ion-molecular chemistry, hydrogenation of simple species on ices, photodissociation and photoionization, freeze-out and desorption, X-ray and cosmic-ray induced reactions, and reactions with vibrationally excited $\rm H_{2}$. The details of these processes are elaborated in Bruderer et al. (2012) and Visser et al. (2018). We used binding energies to grains of 1600 K for CN and 2050 K for HNC and HCN, based on the values from the Kinetic Database for Astrochemistry (KIDA; Wakelam et al. 2012). The chemistry was run in steady-state mode with initial ISM abundances for CNO as listed in Table 1, thus with the ISM C/O. Given the uncertainties in individual rate coefficients and physical structure, our model abundances and column densities have uncertainties of at least a factor of a few. The trends in abundances and line fluxes with model parameters are more robust.

Figure 2 shows the gas densities and abundances (with respect to the total gas density) for HNC and HCN in our fiducial T Tauri disk model. In the hot inner disk ($<$0.5 au), where CN molecules would be all converted into HCN through reaction with $\rm H_{2}$ (Walsh et al., 2015), HCN is predicted to be orders of magnitude more abundant than the range covered in Figure 2. We focus on the outer disk where line emission at (sub-)millimeter wavelengths is generated mostly.

Our model suggests similar spatial distributions for HNC and HCN, with regions of higher HNC and HCN abundance and number density located near the disk surface and extending in the midplane from $\sim$200 au outward (see also Figure 4 in Visser et al. 2018). Their distributions are expected to closely follow the joint distribution of HCNH⁺ and electrons (see Figure 2) because they are mainly produced through dissociative recombination of HCNH⁺. At the surface of the disk, HCN and HNC are photodissociated into H and CN, explaining the upper boundaries on the HCN and HNC abundance distributions.

Although the HCN and HNC dsitributions are largely similar in our model, there are also important differences. Figure 3 shows the 2D HNC-to-HCN abundance ratio distribution for regions in which the HNC gas number density exceeds $\rm 5\times 10^{-5}\,cm^{-3}$. Our model predicts overall low abundance ratios ($<0.3-0.4$) across the disk, with localized higher ratios seen in the cold outer midplane and the warm intermediate layer around 300 au. In the warm layer, the abundance ratio decreases inward toward the hotter region. This temperature-dependent abundance ratio is expected as the destruction of HNC through reaction 4 is more efficient in the warmer regions, leading to lower HNC-to-HCN ratios. The abundance difference between HNC and HCN is also well observed in the column density profiles (see Figure 4), in which the HCN profile is predicted to peak in the inner disk, while the HNC profile shows a double-peak structure.

For a given temperature structure and abundance profile, spectral line emission is obtained through ray-tracing. In this work, we focus on rotational transitions of $J=1-0$ and $J=3-2$ for HNC and HCN, which are close to 90 and 270 GHz, respectively. These lines are accessible with ALMA at the 1.3 and 3 mm bands, within which most disk observations are conducted. Figure 4 shows the velocity-integrated intensity maps for HNC and HCN at the two transitions from our fiducial model, convolved with a $0\aas@@fstack{\prime\prime}2$ Gaussian beam. The disk is set to be face-on (inclination of 0$\degr$) and placed at a distance of 150 pc.

The simulated line emission from both molecules shows a ring-like structure with an emission deficit in the inner disk in the $1-0$ and $3-2$ transitions. Both molecules present a double-ring in both transitions. The relative intensities of the inner and outer rings vary dramatically between species and transitions, however. For HNC lines, emission peaks at 100 au for the $3-2$ transition and at $\sim$300 au for $1-0$ transition; each emission peak corresponds to one of the two peaks in column density profile of HNC. An inner ring is visible for the $1-0$ line and an outer ring for the $3-2$ line, but these are faint compared to the radial emission peaks. For HCN, our model predicts a prominent ring around 100 au for both transitions, while a fainter second ring around 300 au also emerges in the $1-0$ transition. The $3-2$ line presents a small shoulder at 300 au.

The low-amplitude higher-$J$ emission in the outer disk is likely due to high critical densities of their corresponding transitions (see also the discussion of CN emission in Cazzoletti et al. 2018). The critical density of the HNC $3-2$ transition ($\rm\sim 8\times 10^{6}\,cm^{-3}$ at 20 K, typical gas temperature around 300 au) is a few times higher than the local gas density in the outer disk of our fiducial model (see Figure 1). We therefore expect the upper level of this transition to be depopulated, resulting in weak line emission (see Figure 11 for the line emission regions). The lower $J$ transitions with lower critical density ($\rm\sim 3\times 10^{5}\,cm^{-3}$) are still well populated in the low-density outer disk, however, leading to prominent emission there. In addition, the $3-2$ line has a higher upper energy level, which should be substantially less excited in the outer colder disk regions.

The predicted ring-like emission patterns, seen clearly in radial intensity profiles, are consistent with the column density profiles in which peaks of column densities are located off-center (see Figure 4, right panels). Ring-like emission for HNC and HCN is naturally produced in our full-disk models (with a centrally peaked surface density profile) because they cannot be abundantly formed in the dense inner midplane region where the electron density is low (see Figure 2). Molecular rings thus do not necessarily correspond to dust rings that are frequently seen in continuum observations. The ring-like structures of HNC and HCN emission seen in the fiducial model are also present in models with smaller characteristic radius ($R_{\rm c}$ = 15, 30 au), but the ring peak is shifted inward (see Figure 12 for intensity profiles from models with different $R_{\rm c}$). This can be understood when we consider the disk environments that produce the HNC and HCN outer rings: low-density gas where both HCNH⁺ and electrons are present.

Figure 5 shows the predicted radial profiles of HNC-to-HCN line ratio using the intensity profile cuts from the model images. The HNC-to-HCN line ratio generally increases outward, varying from 0.1 to 0.6 in regions where most line emission originates (50 to 350 au in the fiducial model). This increasing pattern is consistently seen in our models with different disk masses (Figure 5). The small bump around 200 au in the 0.1 M_⊙ disk model might arise because disks with higher disk masses are in general cooler, which slows the HNC destruction down.

The radial profiles of the HNC-to-HCN line ratio from our models are consistent with a scenario in which the line ratio is regulated by temperature in disks because the disk temperature decreases outward. To explore this scenario, we investigated the relation between disk temperature and the HNC-to-HCN ratio more directly. Figure 6 shows the extracted HNC-to-HCN abundance ratio and gas temperature at each position of the disk for all disk models with the mass of 0.01 M_⊙. We focus on disk radii outside of 10 au where most of HNC and HCN millimeter emission originates. About half of the parameter space is not occupied; high HNC-to-HCN abundance ratios are not expected in the warmest disk regions. Furthermore, the highest achieved HNC-to-HCN abundance ratio in each temperature bin is predicted to decrease with increasing gas temperature. There is no one-to-one relation between temperature and HNC-to-HCN ratio, however. Figure 3 shows that low ratios ($<0.2$) are expected to appear across the disk, corresponding to a wide range of gas temperatures, and high ratios ($>0.4$) would emerge in the outer disk from regions close to the cold midplane and from the warm layer near the disk surface (note the z/r encoded by the color scheme in Figure 6). We discuss the possible origin for this distribution in Section 4.

Integrated line fluxes are the most readily obtained observables. Thus it is instructive to explore the dependences of line fluxes on key disk properties. Figure 7 shows the predicted general increasing trend of the integrated HNC line fluxes with higher total disk masses. Disks around Herbig stars, with higher UV radiation, are predicted to have brighter (high-$J$) line emission than those around T Tauri stars. The reaction between atomic N and vibrationally excited H₂ (via FUV pumping), initializing one major pathway for the HCNH⁺ formation (Visser et al., 2018), has a strong dependence on UV fluxes, which could then affect the production of HNC and HCN. Our results presented here focus on HNC lines because similar patterns are also expected for HCN lines.

A disk with a larger flaring angle will have more extended disk surface areas exposed to stellar radiation and therefore to UV radiation, which should boost the cyanide chemistry. Figure 7 shows that increasing the disk flaring angle results in higher line fluxes of HNC $3-2$. Our model also predicts that the effect of disk flaring (UV radiation) on the emission of HNC $1-0$ line should be negligible, which can be understood when we consider that most of the $1-0$ emission originates closer to the cold disk midplane in the outer disk, which is barely affected by stellar UV radiation (see Figure 11 in the Appendix).

The UV fluxes in our models have contributions from the stellar blackbody spectrum and from an excess UV emission due to accretion. In addition to the default accretion rate of $10^{-8}\,\rm M_{\odot}\,yr^{-1}$, we ran a set of additional models with different accretion rates ($10^{-7},10^{-9},$ and $10^{-10}\,\rm M_{\odot}\,yr^{-1}$, with order-of-magnitude differences in FUV fluxes, see Visser et al. 2018) to explore the effect of UV fluxes on HNC lines. Similar to what is seen for changing the disk flaring angles, the simulated line fluxes of the HNC $3-2$ transition increase with higher accretion rates, while the fluxes of the HNC $1-0$ lines stay constant (Figure 13). The model-predicted inner ring location for $J=3-2$ line moves outward with increasing UV radiation, consistent with the behavior of CN $N=3-2$ lines as discussed in Cazzoletti et al. (2018), in which the optimal balance between UV fluxes and gas density is reached farther out when UV radiation is stronger. This suggests that the disk UV environment could be probed by observations of higher-$J$ HNC and HCN line transitions.

Recent observations have suggested that substantial fractions of the volatile carbon and oxygen are missing in disks, as traced by CO in the submillimeter (e.g., Favre et al., 2013; Miotello et al., 2017; Long et al., 2017; Zhang et al., 2019) and H₂O vapor in the FIR (Hogerheijde et al., 2011; Du et al., 2017), making the gas-phase carbon and oxygen abundances largely uncertain. In addition, in contrast to an ISM-like C/O ($\sim$0.4, similar to the setup in our fiducial model), the bright C₂H emission in disks requires further oxygen depletion with a C/O$>$1 (Bergin et al., 2016; Miotello et al., 2019). To explore the HNC and HCN line emission under these conditions, we performed three additional carbon- and oxygen-depleted models: (1) The initial overall abundances of carbon and oxygen were depleted by a factor of 10; (2) abundances of carbon and oxygen were depleted by a factor of 100; and (3) under the condition of model (2), oxygen was further depleted to reach a C/O=1.5, all with respect to the fiducial model. The nitrogen abundance was kept the same as in the original model.

Figure 8 presents the column density profiles of HNC and HCN in the depletion models compared to our fiducial model. Our model indicates higher column densities with an increasing level of C and O depletion, especially in the inner disks, and the corresponding line emission is expected to be brighter by a factor of a few. We suspect that the removal of CO modifies the ionization structure in the disk, which then favors the HCN and HNC formation route starting from N₂ and He⁺ (C.12 in Visser et al. 2018). Meanwhile, the destruction of CN with O would slow down with less initial oxygen, leading to a longer CN lifetime and an increased production of HCN and HNC through CN + H₂ reaction. The increase in HCN column densities with enhanced C/O as predicted by our models is consistent with the chemical model results from Cleeves et al. (2018). With higher C/O, our model also predicts more centrally peaked column density profile. The observed diverse morphologies of HCN emission (ring-like or centrally peaked, Bergner et al. 2019) might be explained by different C/O in individual disks.

The HNC-to-HCN line ratio has been shown to correlate with temperature in the ISM. The observed line ratio in cold dark clouds is much higher than what is found in the warmer giant molecular clouds (Irvine & Schloerb, 1984; Goldsmith et al., 1986). A systematic study of the HNC and HCN distribution toward the OMC-1 region found that the HNC-to-HCN ratio is far below 0.1 around Orion-KL, but can reach unity in the colder cloud ridge (Schilke et al., 1992). Recently, Hacar et al. (2020) carried out a dedicated investigation of the HNC-to-HCN ratio in the integral shape filament of Orion and established a strong correlation between the HNC-to-HCN line ratio and the gas kinetic temperature (see the gray curves in Figure 6). The HNC-to-HCN ratio has therefore been proposed to have a great potential as a thermometer in interstellar and circumstellar environments. The predicted increasing trend of HNC-to-HCN line ratio with disk radius suggests that the HNC-to-HCN ratio can be used as probe of disk temperature. However, as discussed in Section 3.2 and Figure 6, a one-to-one correspondence between disk gas temperature and HNC-to-HCN ratio is more challenging to establish.

The HNC-to-HCN abundance ratio map shown in Figure 3 presents a vertical double-peak pattern that cannot simply be explained by temperature-dependent HNC destruction because the gas temperature increases monotonically toward the surface along the vertical direction. To identify other possible explanations, we examined the abundance vertical cuts at 300 au for species involved in HNC/HCN formation and destruction (Figure 10). The disk is assumed to be cold ($\sim$10–20 K) up to 20 au above the disk midplane. In this disk region, the reaction rate for the O destruction pathway should rise with increasing temperature, producing a decreasing HNC-to-HCN ratio when moving upward. Above 20 au, our model predicts a sharp drop in HNC abundance, followed by a gradual increase upward of 25 au. The sudden decrease in HNC abundance (and HNC/HCN) could be explained by the steep increase in O atom abundance around 20 au, while we suspect the following increase in HNC abundace (and HNC/HCN) above a height of 25 au may relate to photodissociation, which gradually takes over as the main destruction pathway for HCN and HNC. Because similar photodissociation cross sections are assumed in the model for the two molecules, we expect a HNC-to-HCN ratio that gradually approaches unity in photodissociation-dominated regions. The CO abundance profile, which starts to decrease above 50 au due to photodissociation, provides the supporting evidence for this hypothesis because the photodissociation of HNC and HCN should occur deeper in the disk than CO photodissociation.

Because the destruction of HNC largely involves C and O, the HNC-to-HCN ratios would be affected by the volatile C and O abundances. In the depletion models discussed in Section 3.4, the consumption of HNC through reactions with C and O is suppressed, and the HNC-to-HCN abundance ratio increases with higher level of elemental depletion in regions close to disk midplane (see Figure 15 in the Appendix, and also the right panel of Figure 8). Meanwhile, high HNC-to-HCN abundance ratios in the intermediate disk layer gradually vanish because the gas temperature increases in the depletion models more rapidly (Figure 10) and the activation of HNC reaction with H would begin deeper in the lower disk layers.

In summary, although the HNC-to-HCN abundance ratio depends on gas temperature in disks, using this ratio as a thermometer in disks will be difficult because the ratio also strongly depends on other parameters, including UV penetration and initial C and O abundances. However, observed abundance ratios above $\sim$0.4 would very likely probe cold or lukewarm ($<$50 K) disk material. A promising way to isolate higher HNC-to-HCN ratios due to low temperatures rather than photodissociation-dominated chemistry is to observe disks with moderate to high inclinations for which the midplane and elevated emitting region could be separated spatially and spectrally. With such observations, the radially resolved HNC-to-HCN ratios in the disk midplane could be deployed as a thermometer, similar to its use in molecular clouds.

Currently available measurements in TW Hya and HD 163296 return disk-averaged HNC-to-HCN $J=3-2$ line ratios of 0.1–0.2 (Graninger et al., 2015). Emission from these shallow observations is therefore likely dominated by the prominent HNC and HCN $J=3-2$ emission component at smaller disk radii. Deep high-resolution observations are required to map the disk radial and vertical thermal structures. The disk around HD 163296, with an intermediate inclination angle, would be an ideal target for this exploration with spatially resolved HCN and HNC observations.

The determination of physical properties of the emitting gas depends on the excitation condition of the lines employed. When the emission lines are thermalized and kinetic temperature is known, molecular column density can be derived from the optically thin lines. If the LTE condition is not satisfied, the inference of physical quantities would not be so straightforward and detailed radiative transfer calculation with model assumptions are required (van Zadelhoff et al., 2001; Teague et al., 2018). As demonstrated in Section 3.2, some of the $3-2$ and $1-0$ emission arises from regions in the outer disk, where the critical densities of HNC and HCN are expected to be higher than the local gas density. We find that the excitation temperature for the $1-0$ line in our model is $\sim 50\%$ of the gas temperature. While the $1-0$ line may be slightly subthermal, the $3-2$ line at this emitting region would be substantially subthermal because this transition has an even higher critical density. Analysis with HNC and HCN lines should therefore take the non-LTE excitation effects into account for more precise constraints of the disk physical conditions.