© ESO, 2018

1. Introduction

Rotating disks of dust and gas around young stars are the cradles of planets, but detailed studies of these objects have only been possible in the last few years since the advent of the Atacama Large Millimeter/submillimeter Array (ALMA). The resolution and sensitivity of the data in the pre-ALMA era only allowed a handful of disks to be characterized, most of the time only in dust (Williams & Cieza 2011; Andrews 2015). The ALMA radio telescope now opens up the possibility to survey hundreds of disks and spatially resolve these disks in gas and dust (Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016).

One particularly intriguing feature of the new images are the ring-like structures observed in dust (Muto et al. 2012; Andrews et al. 2016; Isella et al. 2016) and gas (Bruderer et al. 2014; Kastner et al. 2015; Öberg et al. 2015; Bergin et al. 2016). In some cases, such as the so-called transitional disks (TDs) with inner cavities, a depletion of material in the central regions of the disk can lead to ring-shaped structures for both solids and gaseous molecules (e.g., van der Marel et al. 2016). However, it is apparent that gas and dust can have different distributions and that molecular rings can be due to peculiar physical and chemical conditions rather than to the morphological structure of the disk. In this sense, the morphological structure of the disk and intensity of molecular line emission can provide a wealth of information that cannot be obtained by dust emission alone.

In this scenario, CN is a particularly interesting molecule. One of the brightest molecules after ¹²CO, with fluxes comparable to those of ¹³CO, CN has been observed and studied in disks both with single-dish observations (e.g. Thi et al. 2004; Kastner et al. 2008; Salter et al. 2011; Kastner et al. 2014) and interferometric data (e.g., Dutrey et al. 1997; Öberg et al. 2010, 2011; Chapillon et al. 2012; Guilloteau et al. 2014; Teague et al. 2016; van Terwisga et al. 2017). Although none of the pre-ALMA era data had enough spatial resolution to give information about the morphology of CN emission, CN has been used as a dynamical mass tracer (Guilloteau et al. 2014).

Fig. 1 CN (N = 2−1) ring observed at 0.9′′ (47 AU) in TW Hya with ALMA (data from Teague et al. 2016). The white ellipse at the bottom-left corner shows the synthesized beam size.

CN emission is expected to be sensitive to UV radiation, since its formation relies on the existence of atomic C and N or can result from photodissociation of HCN. Modelling of CN in photodissociation regions (PDRs; e.g. Jansen et al. 1995; Sternberg & Dalgarno 1995; Walsh et al. 2010) and disks (e.g. van Zadelhoff et al. 2003) has been performed and confirms this expectation, placing CN in the upper warm molecular layers. For these reasons, CN is also referred to as a PDR tracer and could thus be a sensitive probe of disk structure.

With ALMA, more information about this molecule is rapidly gathered. The CN (N = 2−1) and CN (N = 3−2) transitions are observable with Band 6 and Band 7, respectively, and the emission can now be spatially resolved. Observations at both wavelengths have provided new information, in particular about the morphology. Teague et al. (2016) used CN (N = 2−1) hyperfine transitions data to measure the velocity dispersion in the TW Hya disk. At the same time, the beam size of 0.5′′ × 0.42′′ allowed Teague et al. (2016) to clearly resolve and distinguish a ring-like structure in their data (Fig. 1). CN was also detected in many disks the ALMA Lupus disk survey (Ansdell et al. 2016). Two of these disks, namely Sz 98 and Sz 71, were large and bright enough to show resolved ring structures, as presented in van Terwisga et al. (2017).

Ring-shaped emission has also been observed for other PDR tracers, such as C₂H. Bergin et al. (2016) concluded that in order to have ring-shaped C₂H emission in full (i.e., non-transitional) disks, a [C]/[O] > 1 ratio is required, which at the same time also explains the high observed fluxes (Kama et al. 2016). In this work, we investigate whether ring-shaped emission is an intrinsic feature of CN in all protoplanetary disks or whether it is due to the specific underlying morphology of some particular disks, such as transitional disks. By modelling CN chemistry we also want to gain further insight into what kind of physical and chemical information about the disk structure can be inferred from resolved CN observations, i.e. to which parameters the CN flux and ring size are most sensitive.

Section 2 describes the code, adopted model, physical structure, and chemical network. In Sect. 3 the results of our models are presented: the main observed trends and the most important reactions responsible for the ring-shaped CN emission. Finally, in Sect. 4 our results are compared with the data available to date, and the main information that can be obtained from CN is highlighted.

2. Model

In this work, we modelled the disks and CN emission using the 2D thermo-chemical code DALI, which combines radiative transfer, chemistry, thermal balance and ray-tracing calculations. This code has been described in detail and extensively tested with benchmark problems (Bruderer et al. 2012; Bruderer 2013) and its results have been compared to observations in a series of other papers (Bruderer et al. 2012; Fedele et al. 2013; Bruderer et al. 2014; Miotello et al. 2017; Facchini et al. 2017). Given an input gas and dust density structure, the dust radiative transfer problem is solved via a 2D Monte Carlo method to obtain the dust temperature structure T_dust and the mean continuum intensities at wavelengths ranging from the UV to mm. Using an initial guess for the gas temperature, the abundances of the molecular and main atomic species are then calculated with either a time-dependent or steady-state chemical network simulation, and their excitation is computed through non-LTE calculations. Finally, the gas temperature structure T_gas can be obtained from the balance between the cooling and heating processes: the new T_gas is then used to re-calculate the chemistry, and the last two steps are iterated until a self-consistent solution is obtained. Both continuum emission and spectral image cubes are finally created using a ray tracer.

2.1. Adopted physical model

For our disk models we assumed a simple parametrized surface density distribution as in Andrews et al. (2011), following an exponentially tapered power law: ${\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{c}}(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{c}}}{)}^{\mathrm{-}\mathit{\gamma }}\exp [\mathrm{-}(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{c}}}{)}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}]\mathit{,}$(1)where R_c is a characteristic radial scale for the surface density profile and Σ_c sets the normalization of the density profile. This surface density profile corresponds to the so-called self-similar solution to the viscous accretion disk model (Lynden-Bell & Pringle 1974; Hartmann et al. 1998), where the shear viscosity parameter is assumed to depend on the radius via ν ∝ R^{− γ} and not to change in time.

The vertical gas density distribution follows a Gaussian, with an aspect ratio that depends on the radius as $\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{h}}_{\mathrm{c}}(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{c}}}{)}^{\mathit{\psi }}\mathrm{·}$(2)In order to mimic dust settling, two populations of grains are considered following the approach of D’Alessio et al. (2006): a small-grains population (0.005−1 μm) and a large-grains population (1−1000 μm). The small-grains population has a scale height of h, while the scale height of the more settled large-grains population is reduced to χh, where 0 <χ< 1. The fraction of dust surface density distributed in the two populations is controlled by the f_large parameter via the relations Σ_dust,large = f_largeΣ_dust and Σ_dust,small = (1−f_large)Σ_dust. We assumed that most of the dust mass in concentrated into the more settled dust grains (e.g., Testi et al. 2003; Rodmann et al. 2006; Lommen et al. 2009; Ricci et al. 2010). After some preliminary testing, we concluded that varying the f_large settling parameter between 0.9 and 0.99 does not produce any significant difference in CN emission. Consequently, we kept f_large fixed in our grid and assumed f_large = 0.99 in all of our models. Also, negligible variation was observed in both the abundance and emission profiles when the settling parameter χ is varied between 0.1 and 0.4. In this work, χ = 0.2 was used.

Other important parameters that affect the chemistry in the models are the stellar and interstellar radiation fields, and the cosmic-ray ionization rate. For the chemistry of CN, the FUV part of the spectrum (6−13.6 eV) becomes particularly important, since photons in this energy range can be responsible for both the formation of CN (through HCN photodissociation) and for its destruction (e.g., Jansen et al. 1995; Sternberg & Dalgarno 1995; Visser et al. 2017; Heays et al. 2017). We modelled different stellar spectra to study the effect of different UV fluxes by assuming the stars to emit as a black body at a given effective temperature T_eff with an excess UV due to accretion. The stellar X-ray spectrum was modelled as a black body at 7 × 10⁷ K between 1 and 100 keV. The interstellar radiation field and the cosmic microwave background were also accounted for and considered to be isotropic. We set the cosmic-ray ionization rate per H₂ to 5 × 10^-17 s^-1.

The computation was carried out on a spatial grid of 265 cells in the radial direction and 120 in the vertical direction. In the radial direction, 35 cells were logarithmically spaced in the first 20 AU and 230 are linearly spaced at ~2.5 AU resolution in the outer disk. The cells in the vertical direction were linearly spaced and their size scales with the aspect ratio h at a given radius. A spectral grid was also used for the radiative transfer, consisting of 58 wavelength bins extending from 912 Å to 3 mm. The abundance for CN and the other chemical species were then derived by assuming steady-state chemistry for the chemical calculations. Such an approach is justified by the fact that the chemistry in the upper layers is fast, therefore making the relevant timescales much shorter than the average age of a Class II disk.

2.2. Grid of models

The goal of our work is to investigate whether or not the ring-like CN emission observed in spatially resolved disks is a common feature in all protoplanetary disks, or if it is due to the specific global morphology of some particular disks such as TDs. The origin of such rings is also investigated. In addition, it is also clear observationally that CN emission can be very bright and it is important to understand which disk and stellar properties most strongly affect the CN line emission intensity. We ran a small grid of models varying some of the parameters describing the physical structure of the disk, dust distribution, and the stellar spectrum. A summary of the parameters used in the models can be found in Table 1.

Fig. 2 Panel a: abundance of CN in a 10-2M⊙, ψ = 0.3, hc = 0.1 disk surrounding a T Tauri star. The location of the peak value is indicated with the green dot. The shape of the distribution of the outer high-CN-abundance region is representative for all the other models. The white line shows the τ = 0.5 contour of the CN (N = 3−2) emission line. The red line shows the dust τ = 1 layer at the same wavelength. The green dotted line indicates the location of the vertical CN abundance maximum for each radius. Panel b: CN (solid blue line) and HCN (cyan line) column densities calculated for the same model, plotted on a linear rather than logarithmic scale. The HCN column density was multiplied by 2 and is more radially concentrated than CN.

2.3. Chemical network

Our models use the chemical network from Bruderer et al. (2012), expanded and updated by Visser et al. (2017) to produce accurate abundances for HCN, HNC, and CN. In particular, this network includes all reactions from the cyanide chemistry review by Loison et al. (2014) and the updated CN photodissociation rate from Heays et al. (2017). Isotopes are not considered.

The network contains standard gas-phase reactions, photoionization and photodissociation, reactions induced by X-rays and cosmic rays. some reactions with PAHs, and a limited grain-surface chemistry. Full details can be found in Appendix A.3.1 from Bruderer et al. (2012) and Visser et al. (2017). Gas-grain interactions include freeze-out, sublimation, and photodesorption. The binding energy of CN is set to 1600 K, as recommended by the Kinetic Database for Astrochemistry (KIDA; Wakelam et al. 2012).

In the surface layers and outer disk, the UV radiation field is strong enough to pump H₂ into vibrationally excited states, denoted as H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$. This vibrational energy can enable neutral-neutral reactions that would otherwise have insurmountable activation barriers (London 1978; Tielens & Hollenbach 1985; Sternberg & Dalgarno 1995; Bruderer et al. 2012). Of particular importance to our models is the reaction of atomic N with H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ to form NH, which in turn reacts with C or C⁺ to form CN or CN⁺ (see below). The initial step has an activation energy of 12 650 K (Davidson & Hanson 1990) and would normally only proceed at gas temperatures above ~1200 K. However, the UV pumping of H₂ in the surface layers allows the N + H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ reaction to proceed regardless of temperature. As such, it is a crucial first step in the formation of CN. The H₂ pumping rate is taken to be eight times its photodissociation rate (Sternberg et al. 2014).

The chemistry is run in steady-state mode, assuming initial molecular abundances as in Cleeves et al. (2015) but with [C]/[H] increased to 10^-4 from 10^-6 unless specified otherwise. Because the surface layer dominates the CN abundance, its chemistry and emission are not found to be sensitive to the specific period over which the chemistry is run and results are similar in steady state.

3. Results

3.1. CN abundance

CN is often referred to as a PDR tracer, and many studies have shown that in protoplanetary disks this molecule forms and survives mostly in the surface layers exposed to intense UV radiation. The same feature is found in our models, with the highest CN abundance region close to the surface and outer edge of the disks.

Figure 2a shows the CN abundance (with respect to the total gas density) for a typical T Tauri disk. Highest CN abundances are located high up in the disk between 10 and 100 AU and extend down to the midplane at the disk edge at large radii where the external UV also plays a role. Such an abundance distribution in turn translates into a ring of high column densities (Fig. 2b), which is consistent with a ring-shaped emission intensity profile (see Sect. 3.2). In contrast, the HCN abundance and column density peak at small radii (see Fig. 2b).

Figure 3 shows a vertical cut of the abundances of the main species at the radial location of the CN abundance maximum (green dot in Fig. 2a). Again, the CN molecule is located near the surface of the disk and between the hydrogen atomic and molecular hydrogen layer, with HCN located deeper into the disk.

Fig. 3 Vertical cut at the radial location of the CN peak abundance. The CN peak is located between the molecular hydrogen (H2) region and the atomic layer.

Fig. 4 Panel a: radial abundances for the most important reactants in the formation of CN (see reactions (4)and (5)), along the dotted green line in Fig. 2a. Each curve is normalized to 1 at the location of the CN maximum. Panel b: H2 excitation rate, H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ de-excitation rate and their ratio along the dotted green line in Fig. 2a. All quantities are normalized to their value at the CN peak abundance (green dot in Fig. 2a) and refer to the model shown in Fig. 2a. The dashed vertical line indicates the location of the CN peak abundance.

Fig. 5 Ray-traced images of the CN emission rings in models with 10-2M⊙, ψ = 0.3, hc = 0.1 disks around a T Tauri and a Herbig star. The ring in the Herbig disk is clearly more extended than that around the T Tauri star. The images are not convolved with any beam, and the resolution is set by the size of the cells in the spatial gridding of the disk (~2.5 AU). For the ray-tracing, the disks are assumed to be at a distance of 150 pc.

It is worth investigating whether the increase in CN column density towards the outer disk is due to decreasing UV, decreasing temperature, or some combination of parameters. The strongest dependence on temperature is found in the back-reaction $\mathrm{CN}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }\mathrm{HCN}\mathrm{+}\mathrm{H}\mathit{,}$(3)which has an activation temperature of 900 K (Baulch et al. 2005). However, while this reaction is important in increasing HCN in the inner disk, it does not play a role in setting the CN abundance since the newly formed HCN is readily photodissociated back to CN.

The main reactions that dominate the formation of CN are instead $\begin{array}{ccc}\mathrm{N}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{\to }& & \mathrm{NH}\mathrm{+}\mathrm{H}\\ {\mathrm{C}}^{\mathrm{+}}\mathrm{+}\mathrm{NH}\mathrm{\to }& & \mathrm{C}{\mathrm{N}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$CN⁺ proceeds then to CN through a charge transfer reaction with H or through the intermediates HCN⁺ and HCNH⁺ followed by dissociative recombination. Alternatively, NH can react with neutral C directly to form CN. Another less important route to cyanides starts from reactions of N with CH₂ and C₂H. Alternatively, reactions of N⁺ with H₂ lead to NH$\begin{array}{c}\mathrm{+}\\ \mathrm{4}\end{array}$, through a series of ion-molecule reactions, which then reacts with C to form HCNH⁺.

CN is primarily destroyed by photodissociation and by reactions with atomic O. The rate coefficient of the CN + O reaction is taken to have no temperature dependence, as discussed in Visser et al. (2017) and consistent with the KIDA recommendations.

The normalized abundances of the most relevant species are shown in Fig. 4a along the dotted green line shown in Fig. 2a. The peak of the CN abundance corresponds to that of NH. In turn, H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ is required in order for NH to exist. NH also shows a second inner peak at 20 AU. This, however, does not translate into a second CN peak because of the lack of C⁺ in this region, which is required for reaction (5) to happen. The net result is that outside 30 AU the normalized abundance profiles of H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$, NH, and CN are qualitatively similar.

The H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ abundance is controlled by the molecular hydrogen excitation and de-excitation rates shown in Fig. 4b. Moving inward from the CN abundance peak, the de-excitation of H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ increases faster than the excitation due to FUV pumping; this leads to a decrease in x(H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$), and consequently also in x(NH) and x(CN). At larger radii the FUV radiation field is attenuated by the dust, de-excitation turns from collisional to spontaneous emission and starts to dominate again, leading to a decrease in x(H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$) and x(CN). Their ratio (solid red line in Fig. 4b) thus has its maximum right at the location of the CN abundance maximum.

In general, the dependence of CN abundance on the balance between H₂ excitation (via FUV pumping) and de-excitation (via collisions) can be summarized with the local ratio between FUV flux G₀ and gas density n_gas. In fact, FUV pumping is proportional to G₀ × n_gas, while collisional de-excitation scales as ${\mathit{n}}_{\mathrm{gas}}^{\mathrm{2}}$. We observe that the value of G₀/n_gas at the CN abundance peak is constant for all of our models; this is consistent with the UV controlling the abundance of CN. When G₀/n_gas is too low, collisional de-excitation dominates over FUV pumping, thus hindering the formation of H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$. On the other hand, when G₀ is high enough, H₂ gets photodissociated. Photodissociation scales with G₀ × n_gas and formation with ${\mathit{n}}_{\mathrm{gas}}^{\mathrm{2}}$. If G₀/n_gas is too high, photodissociation dominates over formation, thus removing molecular hydrogen and preventing the formation of H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$. The H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ layer is located between the molecular and atomic hydrogen layers, and therefore the maximum of H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ and consequently of CN is located where the correct ratio between G₀ and n_gas is met (see Fig. A.1).

Finally, previous astrochemical models have shown that X-rays can play an important role in the chemistry of CN (Lepp & Dalgarno 1996). However, Stäuber et al. (2005, 2007) find that the effect of X-rays in negligible in environments where a high FUV flux is present. Our results indeed do not show significant differences when the L_X parameter is varied between 10²⁷ and 10³¹ erg s^-1.

Fig. 6 Comparison between the radial CN column density profile (blue) and intensity profile (red) for the two disks shown in Fig. 5. The total gas column density multiplied times 10-10 is also plotted in yellow as a reference.

Fig. 7 Simulated observations of a Mdisk = 10-2M⊙, hc = 0.1, ψ = 0.2 disk around a Herbig star, inclined of 60°. The image was simulated using the tasks simobserve and simanalyze in the CASA software.

Fig. 8 Dependence of the CN integrated flux on the total disk mass for disks around T Tauri (upper panel) and Herbig (lower panel) stars. For each value of ψ, the range highlighted by the shaded area is due to the variation of hc. The fluxes are measured by assuming the disks are at a distance of 150 pc. The value measured for TW Hya and rescaled for the distance of the models is also shown.

Fig. 9 Dependence of the CN ring location on the total disk mass in the models.

3.2. Emission

Fig. 10 Different intensity radial profiles for ψ = 0.1,0.2, and 0.3 in a disk with Mdisk = 10-2M⊙ and hc = 0.1 rotating around a T Tauri star. The intensity of the CN emission increases in more flared disks and moves to slightly larger radii.

Fig. 11 Panel a: integrated flux of the CN (N,J) = (3−2, 7/2−5/2) emission in a T Tauri disk with Mgas = 10-2M⊙, ψ = 0.2, and hc = 0.1 as a function of the accretion rate on the central star. For each accretion rate, the UV luminosity between 91.2 and 110 nm is also shown at the top of the figure. The standard model (Ṁ = 10-8M⊙ yr-1) is marked with a blue square. Panel b: same as panel a, for the CN ring radius.

From the temperature structure and chemical abundances obtained for each model it is possible to ray-trace CN emission. In the CN molecule, each rotational energy level (labelled with the quantum number N) is split into a doublet by the spin-rotation interaction (quantum number J), thus leading to three fine structure transitions. Each of the two sub-levels is in turn split into a triplet by nuclear magnetic interactions (quantum number F), giving rise to hyperfine transitions. In this work we ignored the latter transitions, and discuss only the fine structure transitions. In particular, if not otherwise specified, on the brightest of the three, namely the CN (N,J) = (3−2, 7/2−5/2) transition at 340.25 GHz. Figure 5 shows the ray-traced images for two disks from our grid that have the same physical structure but rotate around a T Tauri and a Herbig star, respectively. In the ray tracing, the disks are set to be face-on to allow the CN emission to be clearly distinguishable. The distance of the disks is assumed to be 150 pc.

Fig. 12 Comparison between (panel a) column density radial profiles and (panel b) radial intensity cuts in models with different Rc values.

Fig. 13 Dependence of (panel a) integrated fluxes and (panel b) ring size in models with different Rc values. The standard model (Rc = 60 AU) is marked with a blue square.

As shown in Fig. 2a, CN is optically thin. Accordingly, for T Tauri stars the radial intensity profile is consistent with the radial column density profiles, and the peak of the column density is located at the same radius as that of the emission peak (Fig. 6). In the disks around Herbig stars, CN forms at larger radii than in T Tauri stars. In addition, in some Herbig models, such as that shown in Figs. 5 and 6, emission can peak at even larger radii than the column density peak. This is because the peak in overall CN column density does not correspond to the peak in column density of the CN molecules in the upper state of the transition we are studying (N = 3, J = 7/2), which is located further out. A detailed explanation of this excitation mechanism can be found in Appendix C.

All other models in our study show similar ring-like structures, thus confirming that ring-like emission is an intrinsic feature of CN emission in full disks. In every case, such a structure is due to the higher abundance of CN at larger radii and is not caused by ring-shaped features in the global disk morphology such as those characteristic in transitional disks.

We also concluded that it is not due to sub-thermal excitation effects, which in previous studies have been shown to potentially play an important role for high dipole moment molecules such as CN (e.g., van Zadelhoff et al. 2003), confirming the conclusions of Hily-Blant et al. (2017) for the TW Hya disk. Nevertheless we find that the high critical density of CN does play a role in enhancing the size of CN rings, especially in Herbig stars.

Finally, given the fact that CN in our models is mostly abundant in the uppermost layers of the disk, CN can be used as a probe for the vertical structure of disks. Figure 7 shows a simulated ALMA observation of CN in an inclined disk around a Herbig star, for five representative channels 0.25 km s^-1 wide: even with a modest resolution of 0.22′′ and with 30 min of on-source integration, the two layers are clearly distingushable.

3.3. Dependence on mass and flaring

We investigated how CN emission is affected by some of the key parameters of our models. First, we studied the dependence of the integrated flux and the location of the ring on the total gas mass of the disks.

Figure 8 shows that both for Herbig and T Tauri stars the CN integrated flux increases as the total disk mass increases. The measured fluxes range from 200 mJy km s^-1 and 9 Jy km s^-1, which is consistent with observed values (e.g. Guilloteau et al. 2013). The dependence of the flux on the mass is steeper at low masses for T Tauri stars and flattens for M_disk> 10^-3M_⊙. For the Herbig stars, the flux generally increases less steeply with mass. For both stellar spectra, more flared disks always show higher intensities.

Figure 9, on the other hand, shows the size of the ring as a function of M_disk. More massive disks have larger rings. The difference between T Tauri and Herbig stars is much more evident here than for the flux, hinting at a strong dependence of the CN abundance on the FUV spectrum. In particular, disks surrounding T Tauri stars show rings with sizes ranging between 30 AU and 100 AU while Herbig disks show rings extending to more than 200 AU. The disk scale height, parametrized through h_c, does not play a strong role, as the scatter due to the variation of h_c is small compared to the ranges of the fluxes and of the ring sizes, which from now on we define as the radial location of the emission peak.

Fig. 14 Integrated flux of the CN (N,J) = (3−2, 7/2−5/2) emission in a disk with Mgas = 10-2M⊙, ψ = 0.2, and hc = 0.1 as a function of the C and O depletion. The blue dot refers to the model where a higher gas-phase carbon abundance than oxygen is assumed. All the fluxes are normalized to the non-carbon-depleted δ = 1 case. The fluxes vary by less than a factor of 2.

Finally, disks with a larger flaring angle, exposed to stronger FUV fluxes, show systematically larger rings for both Herbig and T Tauri stars. Changing the flaring angle of the disks in our models has an impact on the CN abundance and consequently on its emission intensity. Figure 10 shows a model with M_disk = 10^-2M_⊙ and h_c = 0.1 rotating around a T Tauri star as an example. The radial intensity cuts are shown for the three different values of ψ used in our grid. As ψ increases, the emission increases in the ring, which moves outward. This confirms our earlier conclusion that UV flux plays an important role in regulating the CN abundance. For each model, the vertical location of the CN peak is determined by the UV field relative to the gas density (Fig. A.1): when the flaring angle is larger the optimal ratio between UV flux and gas density is found at lower altitudes where the gas density is higher.

The same qualitative trends are observed in our models for the (N,J) = (3−2,5/2−5/2) and (N,J) = (3−2,5/2−3/2) fine transitions and for the brightest of the N = 2−1 transitions, namely the J = 5/3−3/2 one at 226.88 GHz. Finally, the ratio between the integrated flux of the (N,J) = (3−2, 7/2−5/2) transition and that of the (N,J) = (2−1,5/2−3/2) one lies between 0.5 and 3 in all our models.

Finally the dependence of the CN emission pattern on the settling and degree of growth of the dust grains was studied. For this purpose, two additional models were run of a 10^-2M_⊙, ψ = 0.2, h_c = 0.1 disk surrounding a T Tauri star assuming f_large = 0.1 and 0.8. In these models, the higher amount of small grains in the disk atmosphere lowers the intensity of the FUV field. Accordingly, the high CN-abundance region is moved to even higher altitudes in the disks in which G₀/n_gas is similar. As a consequence, the CN column density and integrated flux decreases by a factor of ~3 for the f_large = 0.8 model, and the rings are only a few AUs smaller than in the f_large = 0.99 case. In the extreme f_large = 0.1 model, however, the column density and intensity decrease by almost an order of magnitude.

3.4. Excess UV from accretion

In order to further test the impact of the UV field on the CN flux and distribution, we ran an additional set of models by varying the excess UV due to accretion onto a T Tauri star. So far, the T Tauri stellar spectrum was modelled as a 4000 K black-body spectrum with a UV excess for an accretion rate of 10^-8M_⊙ yr^-1. The excess accretion UV was modelled as a 10 000 K black body on a 1.65 R_⊙ radius star (Kama et al. 2016). The additional accretion rates considered are 10^-7, 10^-9, 10^-10M_⊙ yr^-1 and a non-accreting black-body spectrum, and we investigated the effects of the different UV excess on a disk with M_gas = 10^-2M_⊙, ψ = 0.2, and h_c = 0.1.

The integrated flux and ring size as functions of accretion rate are plotted in Fig. 11. The observed trends presented in Fig. 11 confirm the main differences seen between T Tauri and Herbig stars. Higher accretion rates corresponding to higher UV fluxes show higher CN integrated emission. Moreover, the radius of the CN rings increases in size as the UV field becomes stronger. Of course, for a given accretion rate, a spread in the flux and ring size is expected due to the other parameters as explained above.

3.5. Varying the disk size

As a last test, we investigated the effect of different values of R_c on CN morphology, which ultimately reflects on the disk size. As a base model, we used the M_gas = 10^-2M_⊙, ψ = 0.2, h_c = 0.1 T Tauri disk, and then varied R_c keeping the inner surface density profile equal in all the models (see Fig. B.1). In particular, we tested models with R_c = 15,30,60 and 120 AU. Figure 12a shows that the inner part of the CN column density profiles remains constant in all models: this is consistent with the fact that CN abundance is regulated by the balance between stellar FUV and gas density. The CN column density then starts to decrease following the overall gas density profile. Consequently, the CN column density peak is reached at larger radii for higher values of R_c. The same occurs for the radial intensity profile: the ring radius moves increasingly outward as R_c increases from 15 to 60 AU, and becomes almost constant for larger disk radii (Figs. 12b and 13b). Since a high percentage of the intensity comes from the outer disk, the integrated fluxes increase with R_c (Fig. 13a), and the smallest disk shows an order of magnitude lower intensity than the largest. It should be noted that there is a difference in mass between the R_c = 120 AU disk and for the 15 AU of about a factor of ten. However such a difference alone would only provide a variation of a few in flux (Fig. 8a).

3.6. Varying the carbon abundance

Recent observations of protoplanetary disks have shown that lines from CO and its isotopologues, the molecules usually adopted to measure the disk mass, are much weaker than expected in many disks (e.g., Favre et al. 2013; McClure et al. 2016; Miotello et al. 2017; Pascucci et al. 2016; Long et al. 2017; Trapman et al. 2017). This depletion of CO is inferred after photodissociation and freeze-out have been accounted for. One interpretation is that CO is transformed to other less volatile molecules, which is mimicked in our models with an overall volatile C depletion. The gas-phase carbon abundance therefore has recently become one of the largest uncertainties in protoplanetary disks. To account for this effect, we ran three additional models of the M_gas = 10^-2M_⊙, ψ = 0.2, and h_c = 0.1 T Tauri disk. In the first two, the initial overall abundances of carbon, oxygen, sulphur and nitrogen are depleted by a factor of 10 and 100. In the third model, carbon, nitrogen and sulphur are again depleted by a factor of 100, while oxygen is even more depleted in order to get [C] / [O] = 1.5. This assumption is motivated by Kama et al. (2016) and Bergin et al. (2016), where a ratio [C] / [O] > 1 has been used to fit the C₂H flux in TW Hya. As shown in Fig. 14, the CN flux only varies by a factor of less than 2.5 in the four models in spite of two order of magnitude C and O abundance variations, thus showing a weak dependence of CN flux on the level of C and O abundance.

This happens because as the CN abundance in the surface layers decreases, this abundance increases deeper in the disk, because of the lower amount of O, which is the species that CN is mainly destroyed by. For this reason, CN is not a good tracer of C and O depletion in protoplanetary disks.

4. Discussion

The integrated fluxes in our models are consistent with single-dish and interferometric observations. Guilloteau et al. (2013) measured fluxes ranging from a few hundred mJy km s^-1 up to ~10 Jy km s^-1, which is exactly the range of values observed in the models. The same holds for the SMA interferometric observations in Öberg et al. (2010, 2011). Öberg et al. (2011) stated that the radial profiles show different extents for CN and HCN in different disks. In particular, CN appears more radially extended than HCN in V4046 Sgr and LkCa 15. Using more recent ALMA observations, Guzmán et al. (2015) also concluded that CN is ~2 times more extended than HCN in the disk surrounding the Herbig star MWC 480. This feature is well reproduced in our models (see Fig. 2a).

Fig. 15 Comparison between CN, 13CO and C18O integrated flux dependence on the total disk mass. The 13CO and C18O fluxes are taken from Miotello et al. (2017). The solid lines represent the fit functions of the median for 13CO (red) and C18O (blue).

Finally, TW Hya is so far the only disk in which the CN emission has been observed and spatially resolved. CN shows a ring located close to 47 AU radius. With a mass of 0.056 M_⊙ (Bergin et al. 2013) and a flaring index ψ = 0.3 (Kama et al. 2016), both the measured integrated flux for the CN (N,J) = (3−2, 7/2−5/2) transition of 27.1 Jy km s^-1 (Kastner et al. 2014) and the ring size is only slightly lower/smaller than the ψ = 0.3 models we show in Fig. 8a. Still, those values are perfectly consistent with the fact that TW Hya is a low accretor, with Ṁ ~ 10^-9M_⊙ yr^-1 (Debes et al. 2013). For our models we assumed instead Ṁ = 10^-8M_⊙ yr^-1, and, as shown in Fig. 11a, smaller rings and lower fluxes are expected for lower accretion rates.

Recently, Bergin et al. (2016) have shown that in TW Hya other PDR tracers such as C₂H and c-C₃H₂ also present ring-like structures in their emission morphology. In order to reproduce such emission, Bergin et al. (2016) concluded that a strong UV field combined with [C]/[O] > 1 in the upper disk are needed. These conditions are not required for CN, and its ring shape is a natural consequence of the CN formation mechanism and therefore is expected to be a common feature observable in many full disks. The ring and brightness of the CN rings can however be enhanced by a higher UV flux impinging on the disk.

CN emission is therefore a good diagnostic for those parameters affecting UV flux. Among all, flaring is a particularly interesting case. The flaring index is usually derived by fitting the mid-IR range of the SED emitted by the inner few AU of the disk and this value is then assumed to be uniform out to the outer disk. For the first time, information about flaring at much larger radii can be gathered by studying the CN ring radius and flux when some constraints about the spectral type of the central star are available.

CN flux is however not a good tracer for the mass of the disk and CO isotopologues still provide a more sensitive diagnostic (see Fig. 15). This is because in our models CN is most abundant in the upper layers of the disk, and no emission originates from the midplane where the bulk of the gas mass is. Our model results therefore do not support the results of Chapillon et al. (2012), who suggested the existence of large quantities of gas-phase CN close to the disk midplane, based on the low measured values of the excitation temperature with pre-ALMA data.

5. Conclusions

We modelled the CN abundance and emission in full protoplanetary disks around Herbig and T Tauri stars with the 2D thermochemical code DALI to self-consistently solve for the gas temperature structure and molecular abundance. The mass of the disks was varied between 10^-5 and 10^-1M_⊙. We also varied other physical parameters, such as the flaring of the disk, and tested the effect of disk size, UV excess due to accretion onto the central star, and level of depletion of volatile C and O. The modelling shows the following results:

• 1.

  CN shows ring-shaped emission even in full protoplanetarydisks. This is due to the formation route of CN, which is triggeredby the excitation of H₂ to H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ that is mostly abundant in a ring-shaped region close to the surface layer of the disk. H₂ is excited to H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ mostly through UV pumping. For this reason CN emission is very sensitive to UV flux.

• 2.

  Flaring plays an important role: disks with larger ψ show systematically higher fluxes and larger rings. On the other hand, the UV intensity of the central stellar spectrum is critical itself. Herbig stars generally have stronger UV fields, and are thus expected to present larger disks and higher integrated fluxes than T Tauri stars. As for T Tauri stars, UV excess can in first approximation be quantified through the accretion rate (e.g., Kama et al. 2016), with higher accretion rates producing higher UV excess fluxes and consequently brighter and larger CN rings.

• 3.

  The integrated fluxes of the models increase as the total disk mass increases but CN is still less sensitive than the CO isotopologues to the total gas mass. On the other hand, ring size is a much easier parameter to measure and, especially for Herbig stars, it is strongly dependent on the total gas disk mass, therefore providing a useful additional parameter to put constraints on the mass of disks around Herbig stars.

• 4.

  The disk critical radius R_c only plays a role in the CN morphology when smaller than 60 AU: in this case, smaller disks also show smaller CN rings and lower integrated fluxes. However, disks with R_c larger than 60 AU do not show any qualitative difference from disks with R_c = 60 AU. The disk scale height does not strongly affect the CN emission.

The dependence of CN flux and ring size on the above disk parameters make it a sensitive probe of the disk structure. In particular, when the spectral type of the central star is known CN can be combined with the SED to obtain information about parameters that are usually poorly constrained, such as the flaring of the disk at large radii, and to get an independent estimate of the UV flux.

Acknowledgments

We thank the anonymous referee for his/her constructive comments. We also thank A. Miotello, S. van Terwisga, C. Walsh, A. Sternberg, A. Faure, F. Lique and N. van der Marel for useful discussions. Astrochemistry in Leiden is supported by the European Union A-ERC grant 291141 CHEMPLAN, by The Netherlands Research School for Astronomy (NOVA), and by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize. All the figures were generated with the python-based package matplotlib (Hunter 2007).

References

Appendix A: Constant G₀/n_gas

In all of our models, the value of G₀/n_gas is constant to within a factor of a few. Figure A.1 shows this value for our T Tauri models and the same holds for Herbigs. In Fig. A.2 not only the abundance peak, but at each radius the CN abundance maximum (dotted green line) follows the constant G₀/n_gas profile (white line). The model represents a 10^-2M_⊙, ψ = 0.3, h_c = 0.1 around a T Tauri star.

Fig. A.1 Value for G0/ngas where the CN abundance peaks for all the T Tauri models. The different lines for each colour refer to models with different hc.

Fig. A.2 Overlay of the CN abundance structure in a 10-2M⊙, ψ = 0.3, hc = 0.1 disk surrounding a T Tauri star with a green dotted line indicating the location of the vertical CN abundance maximum for each radius. The white line highlights a surface of constant G0/ngas. The two lines follow each other very closely.

Appendix B: Density profiles for different values of R_c

When testing the effect of the disk size on CN emission, we varied the values of R_c keeping the density profiles at the small radii equal for all the models, as shown in Fig. B.1. This keeps the conditions in the inner disk equal for all the models.

Fig. B.1 Comparison between the radial surface density profiles used when comparing models with different Rc parameters. The dotted line indicates, for each profile, the location of Rc.

Appendix C: Level excitation and ring location

As shown in Fig. 6, and as observed also in other disks models around Herbig stars, the CN (N = 3−2) emission peaks at radii which are larger than the radius of the column density peak. This is an excitation effect that originates from the high critical densities of the CN transitions, combined with the low energies of the transitions.

Two main effects play a role:

• 1.

  The CN column density decreases more slowly toward theouter radii for the Herbig disks than for the T Tauriones (see Fig. 6).

• 2.

  The gas density along the CN high-abundance layer decreases with radius (Fig. C.1a). Gas temperature also decreases, but remains always above the energy of the N = 3−2 transitions (~30 K) and also of the N = 5−4 one (~0 K).

What happens is then the following:

• In the inner region of the disk, the gas density and temperature arehigh enough for the lines to be thermalized: as the column densityincreases, the intensity increases accordingly.

• Moving outward along the CN layer, even though the gas temperature remains high, the gas density decreases: as the peak of the CN column density is reached (~ 120 AU in this models) the gas density has fallen below the critical density of the higher level transitions. Consequently, higher level get depopulated and the population of the lower levels grows as we move out and at lower gas densities (see Fig. C.1c).

• Since CN is optically thin, the emission intensity only depends on the population of the upper layers of the transition. The emission thus follows the increase/decrease of populations (Fig. C.1b).

• Moving even further out, at some radius wthe critical density of the N = 3−2 transition is also reached: as explained above level 3 gets unpopulated, the column density of CN in level N = 3 decreases, 3−2 intensity starts decreasing and the population of levels 2 and 1 increases.

This effect is not evident for the T Tauri stars: as mentioned above, for these stars the CN column density profile decreases more steeply, and dominates over the increase in population of the lower levels.

Fig. C.1 Panel a: gas density and temperature along the CN high abundance layer (dotted-green line in Fig. 2a). Panel c: column densities of CN molecules in the (N = 1, J = 3/2) (dotted), (N = 3, J = 7/2) (solid) and (N = 5, J = 11/2) (dashed) levels. Panel b: intensities of the N = 5−4, N = 3−2, and N = 2−1 transitions, for comparison with panel c. The vertical lines mark the location of the column density peaks of CN molecules in the three levels.

All Tables

All Figures

	Fig. 1 CN (N = 2−1) ring observed at 0.9′′ (47 AU) in TW Hya with ALMA (data from Teague et al. 2016). The white ellipse at the bottom-left corner shows the synthesized beam size.
In the text

Fig. 2 Panel a: abundance of CN in a 10-2M⊙, ψ = 0.3, hc = 0.1 disk surrounding a T Tauri star. The location of the peak value is indicated with the green dot. The shape of the distribution of the outer high-CN-abundance region is representative for all the other models. The white line shows the τ = 0.5 contour of the CN (N = 3−2) emission line. The red line shows the dust τ = 1 layer at the same wavelength. The green dotted line indicates the location of the vertical CN abundance maximum for each radius. Panel b: CN (solid blue line) and HCN (cyan line) column densities calculated for the same model, plotted on a linear rather than logarithmic scale. The HCN column density was multiplied by 2 and is more radially concentrated than CN.

In the text

	Fig. 3 Vertical cut at the radial location of the CN peak abundance. The CN peak is located between the molecular hydrogen (H2) region and the atomic layer.
In the text

Fig. 4 Panel a: radial abundances for the most important reactants in the formation of CN (see reactions (4)and (5)), along the dotted green line in Fig. 2a. Each curve is normalized to 1 at the location of the CN maximum. Panel b: H2 excitation rate, H$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ de-excitation rate and their ratio along the dotted green line in Fig. 2a. All quantities are normalized to their value at the CN peak abundance (green dot in Fig. 2a) and refer to the model shown in Fig. 2a. The dashed vertical line indicates the location of the CN peak abundance.

In the text

Fig. 5 Ray-traced images of the CN emission rings in models with 10-2M⊙, ψ = 0.3, hc = 0.1 disks around a T Tauri and a Herbig star. The ring in the Herbig disk is clearly more extended than that around the T Tauri star. The images are not convolved with any beam, and the resolution is set by the size of the cells in the spatial gridding of the disk (~2.5 AU). For the ray-tracing, the disks are assumed to be at a distance of 150 pc.

In the text

	Fig. 6 Comparison between the radial CN column density profile (blue) and intensity profile (red) for the two disks shown in Fig. 5. The total gas column density multiplied times 10-10 is also plotted in yellow as a reference.
In the text

	Fig. 7 Simulated observations of a Mdisk = 10-2M⊙, hc = 0.1, ψ = 0.2 disk around a Herbig star, inclined of 60°. The image was simulated using the tasks simobserve and simanalyze in the CASA software.
In the text

	Fig. 8 Dependence of the CN integrated flux on the total disk mass for disks around T Tauri (upper panel) and Herbig (lower panel) stars. For each value of ψ, the range highlighted by the shaded area is due to the variation of hc. The fluxes are measured by assuming the disks are at a distance of 150 pc. The value measured for TW Hya and rescaled for the distance of the models is also shown.
In the text

	Fig. 9 Dependence of the CN ring location on the total disk mass in the models.
In the text

	Fig. 10 Different intensity radial profiles for ψ = 0.1,0.2, and 0.3 in a disk with Mdisk = 10-2M⊙ and hc = 0.1 rotating around a T Tauri star. The intensity of the CN emission increases in more flared disks and moves to slightly larger radii.
In the text

Fig. 11 Panel a: integrated flux of the CN (N,J) = (3−2, 7/2−5/2) emission in a T Tauri disk with Mgas = 10-2M⊙, ψ = 0.2, and hc = 0.1 as a function of the accretion rate on the central star. For each accretion rate, the UV luminosity between 91.2 and 110 nm is also shown at the top of the figure. The standard model (Ṁ = 10-8M⊙ yr-1) is marked with a blue square. Panel b: same as panel a, for the CN ring radius.

In the text

	Fig. 12 Comparison between (panel a) column density radial profiles and (panel b) radial intensity cuts in models with different Rc values.
In the text

	Fig. 13 Dependence of (panel a) integrated fluxes and (panel b) ring size in models with different Rc values. The standard model (Rc = 60 AU) is marked with a blue square.
In the text

	Fig. 14 Integrated flux of the CN (N,J) = (3−2, 7/2−5/2) emission in a disk with Mgas = 10-2M⊙, ψ = 0.2, and hc = 0.1 as a function of the C and O depletion. The blue dot refers to the model where a higher gas-phase carbon abundance than oxygen is assumed. All the fluxes are normalized to the non-carbon-depleted δ = 1 case. The fluxes vary by less than a factor of 2.
In the text

	Fig. 15 Comparison between CN, 13CO and C18O integrated flux dependence on the total disk mass. The 13CO and C18O fluxes are taken from Miotello et al. (2017). The solid lines represent the fit functions of the median for 13CO (red) and C18O (blue).
In the text

	Fig. A.1 Value for G0/ngas where the CN abundance peaks for all the T Tauri models. The different lines for each colour refer to models with different hc.
In the text

	Fig. A.2 Overlay of the CN abundance structure in a 10-2M⊙, ψ = 0.3, hc = 0.1 disk surrounding a T Tauri star with a green dotted line indicating the location of the vertical CN abundance maximum for each radius. The white line highlights a surface of constant G0/ngas. The two lines follow each other very closely.
In the text

	Fig. B.1 Comparison between the radial surface density profiles used when comparing models with different Rc parameters. The dotted line indicates, for each profile, the location of Rc.
In the text

Fig. C.1 Panel a: gas density and temperature along the CN high abundance layer (dotted-green line in Fig. 2a). Panel c: column densities of CN molecules in the (N = 1, J = 3/2) (dotted), (N = 3, J = 7/2) (solid) and (N = 5, J = 11/2) (dashed) levels. Panel b: intensities of the N = 5−4, N = 3−2, and N = 2−1 transitions, for comparison with panel c. The vertical lines mark the location of the column density peaks of CN molecules in the three levels.

In the text