A study of singly deuterated cyclopropenylidene
c-C3HD in protostar IRAS 16293-2422

Liton Majumdar, P. Gratier, I. Andron, Valentine Wakelam, E. Caux

► To cite this version:

Liton Majumdar, P. Gratier, I. Andron, Valentine Wakelam, E. Caux. A study of singly deuterated cyclopropenylidene c-C3HD in protostar IRAS 16293-2422. Monthly Notices of the Royal Astronomical Society, 2017, 467 (3), pp.3525-3532. 10.1093/mnras/stx259. hal-01451606

HAL Id: hal-01451606

https://hal.science/hal-01451606v1

Submitted on 24 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A study of singly deuterated cyclopropenylidene *c*-C₃HD in the protostar IRAS 16293–2422

L. Majumdar,<fnref n="1" />,<fnref n="2" />* P. Gratier,<fnref n="1" />* I. Andron,<fnref n="1" /> V. Wakelam<fnref n="1" />* and E. Caux<fnref n="3" />,<fnref n="4" />

Accepted 2017 January 27. Received 2017 January 27; in original form 2016 September 15

## ABSTRACT

Cyclic-C₃HD (*c*-C₃HD) is a singly deuterated isotopologue of *c*-C₃H₂, which is one of the most abundant and widespread molecules in our Galaxy. We observed IRAS 16293–2422 in the 3-mm band with a single frequency set-up using the EMIR heterodyne 3 mm receiver of the IRAM 30-m telescope. We observed seven lines of *c*-C₃HD and three lines of *c*-C₃H₂. Observed abundances are compared with astrochemical simulations using the NAUTILUS gas-grain chemical model. Assuming that the size of the protostellar envelope is 3000 au and assuming the same excitation temperatures for both *c*-C₃H₂ and *c*-C₃HD, we obtain a deuterium fraction of 14±4 per cent.

**Key words:** astrochemistry – methods: statistical – ISM: abundances – ISM: evolution – ISM: molecules.

## 1 INTRODUCTION

The cyclic form of C₃H₂, cyclopropenylidene (*c*-C₃H₂), is one of the most abundant and widespread molecules in our Galaxy (Matthews & Irvine 1985). It is also the first hydrocarbon ring molecule found in space. It has been detected in different astrophysical environments: from diffuse gas to cold dark clouds, giant molecular clouds, photodissociation regions, circumstellar envelopes and planetary nebulae (Thaddeus, Vrtilek & Gottlieb 1985; Cox, Guesten & Henkel 1987; Vrtilek, Gottlieb & Thaddeus 1987; Madden et al. 1989; Lucas & Liszt 2000).

*c*-C₃H₂ shows strong rotational transitions in millimetre wavebands because of its large dipole moment and favourably small partition function, which make this molecule a useful probe of physical conditions where it is found (see Bell et al. 1986 for a discussion). *c*-C₃H₂ is often used as a template to study deuteration along different sources as the rotational transitions of its isotopic species are also very intense. *c*-C₃H₂ can exist in singly (*c*-C₃HD) and doubly (*c*-C₃D₂) deuterated forms. *c*-C₃HD was first detected in the dense molecular cloud TMC-1 by Bell et al. (1986). Gerin et al. (1987) measured a 1:5 ratio (i.e. 20 per cent) of the 2₁,₂ → 1₀,₁ lines of *c*-C₃HD to *c*-C₃H₂ in TMC-1. According to Gerin et al. (1987), this high degree of deuterium fractionation indicates that the emission comes from extremely cold dense cores (*T*<sub>rot</sub> = 4–7 K) and thus *c*-C₃HD samples the coldest part of the interstellar medium (ISM). *c*-C₃D₂ was also detected recently by Spezzano et al. (2013)

towards the starless cores TMC-1C and L1544. They measured the abundance of *c*-C₃D₂ with respect to *c*-C₃H₂ and found it to be 0.4–0.8 per cent in TMC-1C and 1.2–2.1 per cent in L1544. This clearly shows that the abundance of *c*-C₃D₂ is much lower than that of *c*-C₃HD and thus supports the idea that the singly deuterated species forms more easily in the ISM compared to the multiply deuterated species (Brown & Millar 1989; Willacy & Millar 1998). This makes *c*-C₃HD a better candidate than *c*-C₃D₂ to study the chemistry and the physical conditions of cold environments. The main goal of the present paper is to report the detection of *c*-C₃HD in the low-mass protostar IRAS 16293–2422 (hereafter IRAS 16293) and to study deuterium fractionation of *c*-C₃H₂ from both observational and theoretical points of view.

This paper is structured as follows. In Section 2, we give a detailed description of our observations along with methodology for the analysis. In Section 3, we discuss the chemical model, which includes our latest deuterium chemical network with associated spin chemistry. Finally, we discuss our results in the last section.

## 2 OBSERVATIONS AND DATA REDUCTION

### 2.1 Observations

The observations were obtained at the Institut de Radioastronomie Millimétrique (IRAM) 30-m telescope during the period 2015 August 18–23. Overall, the weather conditions were similar to an average summer (a median value of 4–6 mm water vapour). We performed our observations by using the heterodyne 3-mm eight mixer receiver (EMIR) tuned at a frequency of 89.98 GHz in the lower inner sideband. This receiver was followed by a Fourier

<footnote n="1">¹Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="2">²Indian Centre for Space Physics, 43 Chalantika, Garia Station Road, Kolkata, 700084, India</footnote>

<footnote n="3">³Université de Toulouse, UPS-OMP, IRAP, Toulouse, France</footnote>

<footnote n="4">⁴CNRS, IRAP, 9 Av. Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote>* E-mail: liton.icsp@gmail.com (LM); pierre.gratier@u-bordeaux.fr (PG); valentine.wakelam@u-bordeaux.fr (VW)</footnote>
<table><caption>**Table 1.** Observed lines and spectroscopic parameters for *c*-C₃HD and *c*-C₃H₂.</caption><thead><tr><th>Species</th><th>Integrated flux<br>(K km s<sup>-1</sup>)</th><th>V<sub>LSR</sub><br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>Observed<br>frequency (MHz)</th><th>A<sub>ij</sub><br>(s<sup>-1</sup>)</th><th>E<sub>up</sub><br>(K)</th><th>Quantum<br>numbers</th></tr></thead><tbody><tr><td rowspan="7">c-C<sub>3</sub>HD<sup>a</sup></td><td>0.013 ± 0.009</td><td>4.3 ± 0.162</td><td>1.0 ± 0.70</td><td>85643.318</td><td>1.43 × 10<sup>-5</sup></td><td>26.6</td><td>4<sub>3,2</sub> → 4<sub>2,3</sub></td></tr><tr><td>0.028 ± 0.003</td><td>4.3 ± 0.064</td><td>1.3 ± 0.12</td><td>102423.019</td><td>1.53 × 10<sup>-5</sup></td><td>22.3</td><td>4<sub>1,3</sub> → 4<sub>0,4</sub></td></tr><tr><td>0.243 ± 0.024</td><td>4.2 ± 0.016</td><td>1.1 ± 0.02</td><td>104187.126</td><td>3.96 × 10<sup>-5</sup></td><td>10.8</td><td>3<sub>0,3</sub> → 2<sub>1,2</sub></td></tr><tr><td>0.047 ± 0.006</td><td>4.2 ± 0.031</td><td>0.8 ± 0.08</td><td>104799.707</td><td>7.29 × 10<sup>-6</sup></td><td>10.9</td><td>3<sub>1,3</sub> → 2<sub>1,2</sub></td></tr><tr><td>0.032 ± 0.005</td><td>4.5 ± 0.086</td><td>1.4 ± 0.19</td><td>106256.108</td><td>1.69 × 10<sup>-5</sup></td><td>22.5</td><td>4<sub>2,3</sub> → 4<sub>1,4</sub></td></tr><tr><td>0.054 ± 0.006</td><td>4.1 ± 0.025</td><td>0.9 ± 0.07</td><td>106811.09</td><td>7.87 × 10<sup>-6</sup></td><td>10.8</td><td>3<sub>0,3</sub> → 2<sub>0,2</sub></td></tr><tr><td>0.280 ± 0.028</td><td>4.1 ± 0.005</td><td>0.8 ± 0.01</td><td>107423.671</td><td>4.47 × 10<sup>-5</sup></td><td>10.9</td><td>3<sub>1,3</sub> → 2<sub>0,2</sub></td></tr><tr><td rowspan="3">c-C<sub>3</sub>H<sub>2</sub><sup>b</sup></td><td>0.125 ± 0.013</td><td>4.2 ± 0.018</td><td>1.4 ± 0.06</td><td>84727.688</td><td>1.04 × 10<sup>-5</sup></td><td>16.14</td><td>3<sub>2,2</sub> → 3<sub>1,3</sub></td></tr><tr><td>1.372 ± 0.137</td><td>4.2 ± 0.010</td><td>1.4 ± 0.006</td><td>85338.894</td><td>2.32 × 10<sup>-5</sup></td><td>6.45</td><td>2<sub>1,2</sub> → 1<sub>0,1</sub></td></tr><tr><td>0.173 ± 0.017</td><td>4.3 ± 0.028</td><td>1.9 ± 0.05</td><td>85656.431</td><td>1.52 × 10<sup>-5</sup></td><td>29.07</td><td>4<sub>3,2</sub> → 4<sub>2,3</sub></td></tr></tbody></table>

*Notes.* <sup>a</sup>Lines from Bogey et al. (1987).

<sup>b</sup>Lines from Thaddeus et al. (1985).

transform spectrometer in its 195-kHz resolution mode. The observed spectrum was composed of two regions: 84.4–92.3 GHz and 101.6–107.9 GHz.

The observations were made towards the midway point between sources A and B of IRAS 16293 at $\alpha_{2000} = 16^h32^m22^s75$, $\delta_{2000} = -24^\circ28'34''/2$. The A and B components, separated by 5.5 arcsec, are both inside the telescope beam of our observations at all frequencies. All observations were performed using the wobbler switching mode with a period of 2 s and a throw of 90 arcsec, ensuring mostly flat baselines even in summer conditions and observations at low elevation. At the beginning of each observing run, the closest planet, Saturn, was used for focus. Pointing was checked every hour with mostly good pointing corrections (less than a third of the telescope beam size 30 arcsec).

## 2.2 Results

### 2.2.1 *c*-C₃HD and *c*-C₃H₂ line properties

The data were reduced and analysed using the CLASS software from the GILDAS<fnref n="1" /> package. We made Gaussian fits to the detected lines following a local low-order (0 or 1) polynomial baseline subtraction. Table 1 shows the result of these fits for the seven observed lines of *c*-C₃HD and three observed lines of *c*-C₃H₂. In Table 1, uncertainties on integrated flux are the Gaussian fit uncertainties (with an added 10 per cent calibration error), uncertainties on $V_{LSR}$ are the quadratic sum of the Gaussian fit statistical uncertainties and the frequency uncertainties from the spectroscopic catalogue, and uncertainties on FWHM are the Gaussian fit uncertainties only. All the observed lines are single-component features with the mean LSR velocity of ~4.2 km s⁻<fnref n="1" /> and the mean FWHM of ~1 km s⁻<fnref n="1" />. Caux et al. (2011) identified four types of kinematical behaviours for different species based on the different FWHM and $V_{LSR}$ distributions. We found that *c*-C₃HD and *c*-C₃H₂ belong to the type I (i.e. FWHM ≤ 2.5 km s⁻<fnref n="1" />, $V_{LSR}$ ~ 4 km s⁻<fnref n="1" />, $E_{up}$ ~ 0–50 K), which corresponds to species abundant in the cold envelope of IRAS 16293.

### 2.2.2 Radiative transfer modelling for *c*-C₃HD and *c*-C₃H₂

Majumdar et al. (2016b) used a Bayesian model to recover the distribution of parameters that best agree with the observed line intensities for CH₃SH in the same source. Here, we use the same

<table><caption>**Table 2.** Distribution functions of priors for the parameters used in the Bayesian approach. Uniform (min-value, max-value) is a uniform distribution with values going from the minimum value to the maximum value.</caption><thead><tr><th>Parameter</th><th>Distribution</th></tr></thead><tbody><tr><td>log<sub>10</sub>N (cm<sup>-2</sup>)</td><td>Uniform (8, 22)</td></tr><tr><td>log<sub>10</sub>T<sub>ex</sub> (K)</td><td>Uniform (log<sub>10</sub>3, log<sub>10</sub>200)</td></tr><tr><td>V (km s<sup>-1</sup>)</td><td>Uniform (2.8, 4.8)</td></tr><tr><td>ΔV (km s<sup>-1</sup>)</td><td>Uniform (0.25, 10)</td></tr><tr><td>log<sub>10</sub>σ<sub>add</sub> (K)</td><td>Uniform (-3, 1)</td></tr></tbody></table>

model with the exception that not only the integrated intensities but also the full modelled spectra are compared to the observed line intensities.

We use the local thermal equilibrium (LTE) radiative transfer code implemented in the GILDAS–Weeds package (Maret et al. 2011) to model the emission of *c*-C₃HD and *c*-C₃H₂. The input parameters in this radiative transfer model are the species column density, the linewidth, the excitation temperature and the source size, along with an accurate spectroscopic catalogue. The spectroscopic catalogues for *c*-C₃HD and *c*-C₃H₂ were retrieved from the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005). All the observed frequencies along with their Einstein coefficients, upper-level energies and the associated quantum numbers are listed in Table 1. Here, we have chosen to fix the source size to 25 arcsec (~3000 au, a typical size of the protostellar envelope by Caux et al. 2011) because the range of frequencies is small. The likelihood function assumes that the errors are normally distributed, with a noise term consisting of a sum in quadrature of the observed per channel uncertainty and an additional noise term left as a free parameter of the model. The priors are chosen to be uniform and non-informative over the range of variation as defined in Table 2.

The sampling of the posterior distribution function is carried out using the PYTHON implementation EMCEE<fnref n="2" /> (Foreman-Mackey et al. 2013) of the Affine Invariant Ensemble Monte Carlo Markov Chain approach (Goodman & Weare 2010). Sixty walkers are initialized in a small six-dimensional sphere in the centre of the parameter space. The chains are evolved for a burn-in sequence of 3000 steps,

<footnote n="1">¹ https://www.iram.fr/IRAMFR/GILDAS/</footnote>

<footnote n="2">² https://github.com/dfm/emcee</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. The top three panels are for c-C₃H₂ and the bottom seven panels are for c-C₃HD. The red lines show the observed lines attributed to c-C₃HD and c-C₃H₂. The black lines show the distribution of modelled spectra following the posterior distribution of parameters shown in Figs 2 and 3. The thick line denotes the median of the distribution. The dark and light grey regions show the 68 and 95 per cent confidence intervals, respectively. The dotted lines are 1σ, 3σ and 5σ noise levels.</figcaption></figure>

after which the convergence is checked by examining a plot of the running mean of each parameter.

Fig. 1 shows the comparison of observed spectra and modelled spectra. Figs 2 and 3 show the one-dimensional and two-dimensional histograms of the posterior probability distribution function for c-C₃HD and c-C₃H₂, respectively. In both cases, all parameters are well defined.

Table 3 summarizes the one point statistics for the marginalized posterior distributions of parameters. The uncertainties reported in Table 3 are 1σ symmetric error bars. There is about a factor of 2 difference in linewidth between c-C₃HD and c-C₃H₂ (shown in Table 3); this cannot arise from opacity effects as all the lines are optically thin except c-C₃H₂ line at 85.338 GHz. One possibility is that c-C₃H₂ comes from both the central hotter region and the outer envelope, while c-C₃HD only comes from the outer cold envelope. This is further strengthened by the chemical modelling,

which shows such a pattern in Fig. 4. Thus, the deuterium fraction we obtain in Section 2.3 is a lower limit in the external envelope.

## 2.3 Observational constraint on deuterium fraction for c-C₃H₂

The observed main beam temperature $T_{mb}$, for a beam of size $\theta_t$ of a Gaussian source of size $\theta_s$ is

$$ T_{mb} = \eta[J(T_{ex}) - J(T_{CMB})](1 - e^{-\tau}) \\ = \frac{\theta_s^2}{\theta_s^2 + \theta_t^2}[J(T_{ex}) - J(T_{CMB})](1 - e^{-\tau}), \quad (1) $$

where $\eta$ is the beam dilution factor (which is identical for two molecules as we fix the source size), $J(T) = (h\nu/K)1/(e^{h\nu/KT} - 1)$ is the radiation temperature and $\tau$ is the opacity. For optically thin conditions, we can assume that $\tau \propto N$. Under this assumption, $T_{mb} \propto [J(T_{ex}) - J(T_{CMB})]N$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** One-dimensional and two-dimensional histograms of the posterior distribution of parameters for c-C$_3$H$_2$. Contours contain 68 and 95 per cent of samples, respectively. The quoted uncertainties are statistical only, without the 10 per cent calibration error.</figcaption></figure>

Under optically thin conditions, we can assume that the excitation properties of the two isotopologues are going to be similar, which is confirmed by the fact that we find similar values for the excitation temperatures when carrying out an independent analysis. We need more lines, and therefore more sensitive observations, to see whether temperatures are in fact different. The choice to use the same temperature has also been made in the past for the LTE modelling of isotopologues; see, for instance, Caselli et al. (2002) and Daniel et al. (2013).

Now we can represent the column density ($N_2$) of c-C$_3$H$_2$ at the excitation temperature (6.3 K) of c-C$_3$HD as

$$N_2 = N_1[J(T_{ex1}) - J(T_{CMB})]/[J(T_{ex2}) - J(T_{CMB})], \quad (2)$$

where $N_1$ is the column density of c-C$_3$H$_2$ at 8.9 K. At frequency 85 GHz, $J(T_{ex1}) = 7.02$ for $T_{ex1}$ of 8.9 K and $J(T_{ex2}) = 4.46$ for

$T_{ex2}$ of 6.3 K [$J(T_{CMB}) = 1.17$ at $T_{CMB} = 2.73$ K], and we derive a column density of $N(c-C_3H_2) = N_2 = 1.78N_1 = 4.7 \times 10^{13}$ cm$^{-2}$. Thus, we obtain a deuteration fraction $N(c-C_3HD)/N(c-C_3H_2) = 6.7 \times 10^{12}/4.7 \times 10^{13} = 14$ per cent.

If we define the errors on column densities of c-C$_3$H$_2$ and c-C$_3$HD by $errN_2$ and $errN_1$, the statistical errors on column densities of c-C$_3$H$_2$ and c-C$_3$HD by $errstatN_2$ and $errstatN_1$ and relative calibration error by $errcal$ (10 per cent in both cases), then we can write

$$\frac{errN_2}{N_2} = \sqrt{\left(\frac{errstatN_2}{N_2}\right)^2 + (errcal)^2} \\ = \sqrt{0.02^2 + 0.1^2} = 0.10 \quad (3)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** One-dimensional and two-dimensional histograms of the posterior distribution of parameters for c-C₃HD. Contours contain 68 and 95 per cent of samples, respectively. The quoted uncertainties are statistical only, without the 10 per cent calibration error.</figcaption></figure>

and

$$ \frac{\text{err}N_1}{N_1} = \sqrt{\left(\frac{\text{errstat}N_1}{N_1}\right)^2 + (\text{errcal})^2} \\ = \sqrt{0.1^2 + 0.1^2} = 0.14. \qquad (4) $$

For N(c-C₃HD)/N(c-C₃H₂) = R, then

$$ \begin{aligned} R_{\text{max}} &= N(c\text{-C}_3\text{HD})_{\text{max}}/N(c\text{-C}_3\text{H}_2)_{\text{min}} \\ &= 6.7 \times (1 + 0.14) \times 10^{12}/4.7 \times (1 - 0.1) \times 10^{13} \\ &= 18\% \end{aligned} \qquad (5) $$

<table><caption>**Table 3.** Point estimates of the posterior distribution function corresponding to the median and 1σ uncertainty. Quoted uncertainties are statistical only, without the 10 per cent calibration error.</caption><thead><tr><th>Parameter</th><th>c-C<sub>3</sub>HD</th><th>c-C<sub>3</sub>H<sub>2</sub></th></tr></thead><tbody><tr><td>log<sub>10</sub>N (cm<sup>-2</sup>)</td><td>12.83 ± 0.05</td><td>13.42 ± 0.01</td></tr><tr><td>log<sub>10</sub>T<sub>ex</sub> (K)</td><td>0.80 ± 0.05</td><td>0.95 ± 0.03</td></tr><tr><td>log<sub>10</sub>ΔV (km s<sup>-1</sup>)</td><td>-0.09 ± 0.02</td><td>0.11 ± 0.01</td></tr><tr><td>log<sub>10</sub>[X]<sup>a</sup></td><td>-11.3 ± 0.05</td><td>-10.7 ± 0.01</td></tr></tbody></table>

Note. <sup>a</sup>[X] = N(X)/N(H₂).

Values with 95 per cent confidence.
<figure><img src="image_6.png" /><figcaption>Figure 4. Abundance with respect to H₂ for c-C₃H₂ and c-C₃HD predicted by our model as a function of radius. sc-C₃H₂ and sc-C₃HD represent c-C₃H₂ and c-C₃HD on the surface of grains, respectively.</figcaption></figure>

and

$$ 
\begin{aligned}
R_{\min} &= N(c-C_3HD)_{\min}/N(c-C_3H_2)_{\max} \\
&= 6.7 \times (1 - 0.14) \times 10^{12}/4.7 \times (1 + 0.1) \times 10^{13} \\
&= 11\%.
\end{aligned} 
\quad (6) $$

Thus, the error on the deuterium fraction can be written as $N(c-C_3HD)/N(c-C_3H_2) = 14_{-3}^{+4}$ per cent. Bell et al. (1988) also measured deuterium fractionation of c-C₃H₂ in the order of 5–15 per cent in TMC-1.

# 3 MODELLING OF THE PROTOSTELLAR ENVELOPE

## 3.1 NAUTILUS chemical model with deuteration

To model the chemistry of c-C₃HD and c-C₃H₂ in IRAS 16293, we used the same approach as Majumdar et al. (2016b). For our simulation, we have used the two-phase version of the NAUTILUS gas-grain chemical model (Majumdar et al. 2016b; Ruaud et al. 2015; Wakelam et al. 2015) with deuteration and spin chemistry (Majumdar et al., in preparation). NAUTILUS allows the computation of the chemical composition as a function of time in the gas phase and at the surface of interstellar grains. All the equations and the chemical processes included in the model are described in detail in Ruaud et al. (2015). In the current model, several types of chemical reactions are considered in the gas phase by following the kida.uva.2014 chemical network of Wakelam et al. (2015) with the recent extension to deuteration and spin chemistry of Majumdar et al. (2016a). In the two-phase version of NAUTILUS, there is no differentiation between the species in the mantle and at the surface. NAUTILUS considers interaction between gas and grains via four major processes: physisorption of gas-phase species on to grain surfaces, diffusion of the accreted species, reaction at the grain surface and finally evaporation to the gas phase. Our model also considers different types of evaporation processes such as thermal evaporation, evaporation induced by cosmic rays (following Hasegawa & Herbst 1993) and chemical desorption, as suggested by Garrod, Wakelam & Herbst (2007). We adopt similar initial elemental abundances as those reported

in Hincelin et al. (2011) with a deuterium and fluorine elemental abundance relative to hydrogen of $1.6 \times 10^{-5}$ (Linsky et al. 2006) and $6.68 \times 10^{-9}$ (Neufeld, Wolfire & Schilke 2005), respectively. Here, the species are assumed to be initially in an atomic form as in diffuse clouds, except for hydrogen and deuterium, which are initially in H₂ and HD forms, respectively. For our standard model, we have used a C/O ratio of 0.7 and an ortho-to-para H₂ ratio of 3.

In our current model, all the variables relative to H₂ have been modified in terms of nuclear spin states ortho-H₂ and para-H₂. It has been known for a long time that H₂, D₂, H₃⁺, H₂D⁺, D₂H⁺ and D₃⁺, along with their spin isomers, are the main species that dictate deuterium fractionation at low temperatures (see Ceccarelli et al. 2014). A detailed description and benchmarking of our deuterated network with spin chemistry are presented in Majumdar et al. (2016a). This network is available on the KIDA<fnref n="3" /> web site. Our network for surface reactions and gas-grain interactions is based on the one from Garrod et al. (2007) with added deuteration and spin chemistry from Majumdar et al. (2016a).

## 3.2 One-dimensional physical structure

To follow the deuterium fractionation of c-C₃H₂ in IRAS 16293, we have used the same one-dimensional physical structure as in Aikawa et al. (2008), Wakelam et al. (2014) and Majumdar et al. (2016b). This physical structure was based on the non-grey radiation hydrodynamic model by Masunaga & Inutsuka (2000) to follow the core evolution from the pre-stellar core to the protostellar core. It starts from a hydrostatic pre-stellar core with a central density $n(H_2) \sim 3 \times 10^4$ cm⁻<fnref n="3" />. The core is extended up to $r = 4 \times 10^4$ au with a total mass of $3.852 M_{\odot}$, which exceeds the critical mass for gravitational instability.

The initial temperature for the core is around 7 K at the centre and around 8 K at the outer edge. Here, cosmic ray heating, cosmic background radiation and ambient stellar radiation balance the cooling caused by dust thermal emission. In the model, the core

<footnote n="3">³ http://kida obs.u-bordeaux1.fr/</footnote>
stays at its hydrostatic structure for $1 \times 10^6$ yr to set up the initial molecular conditions for the collapse stage. After $1 \times 10^6$ yr, the core starts to contract and the contraction is almost isothermal as long as the cooling is efficient. Eventually, the compressional heating overwhelms the cooling, which causes a rise in temperature in the central region. Contraction then decelerates because of the increase in the gas pressure, which eventually makes the first hydrostatic core, known as the first core at the centre. When the core centre reaches very high density ($10^7$ cm$^{-3}$) and high temperature (2000 K), the hydrostatic core becomes unstable as a result of H$_2$ dissociation and starts to collapse again. This collapse is referred to as the second collapse. Within a short period of time, the dissociation degree approaches unity at the centre because of the rapid increase of central density. Then, the second collapse ceases, and the second hydrostatic core (i.e. the protostar) is formed and the infalling envelope around this protostar is known as the protostellar core. In the model, the initial pre-stellar core evolves to the protostellar core in $2.5 \times 10^5$ yr. When the protostar is formed, the model further follows the evolution for $9.3 \times 10^4$ yr, during which the protostar grows by mass accretion from the envelope.

## 4 MODELLING RESULTS AND DISCUSSION

Here, we present the predicted abundances of c-C$_3$HD and c-C$_3$H$_2$ in the protostellar envelope using the physical and chemical models previously described, and we compare with our observations. Fig. 4 shows the computed abundances of c-C$_3$HD and c-C$_3$H$_2$, in the gas phase and at the surface of the grains in the protostellar envelope as a function of radii, at the end of the simulations (i.e. for a protostellar age of $9.3 \times 10^4$ yr). We considered abundances at $9.3 \times 10^4$ yr because the physical structure of the envelope at this age is similar to the one constrained in the envelope of IRAS16293 by Crimier et al. (2010) from multiwavelength dust and molecular observations (see Wakelam et al. 2014 for a discussion).

The abundance profiles of c-C$_3$HD and c-C$_3$H$_2$ predicted by our model can be divided into three regions. The first region is defined by radii larger than 200 au and temperatures below 50 K. In this region, the gas-phase abundance of both c-C$_3$HD and c-C$_3$H$_2$ decrease towards the centre of the envelope due to high depletion because of density increase. In the outer part of the envelope (radii greater than 2000 au and temperatures below 30 K), gas-phase c-C$_3$HD forms mainly by the dissociative recombination of c-C$_3$H$_2$D$^+$, c-C$_3$H$_2$D$^+$ is produced mainly from the deuteron transfer from para-H$_2$D$^+$ and DCO$^+$ to c-C$_3$H$_2$. c-C$_3$H$_2$D$^+$ is also partly produced by c-C$_3$H$_2$ + H$_2$DO$^+$ and para-H$_2$ + C$_3$D$^+$ reactions. So, in the cold outer envelope, deuteration of c-C$_3$HD is mainly controlled by para-H$_2$D$^+$ and DCO$^+$. In these regions, c-C$_3$H$_2$ also mainly forms from the dissociative recombination reaction of c-C$_3$H$_3^+$. From 2000 to 200 au (where the temperature started to increase from 30 K to 50 K), c-C$_3$HD mainly forms from the neutral-neutral reaction CH + C$_2$HD → H + c-C$_3$HD whereas c-C$_3$H$_2$ mainly forms from the H + CH$_2$CCH → para-H$_2$ + c-C$_3$H$_2$ reaction.

Between 200 and 100 au, the gas-phase abundance of both c-C$_3$HD and c-C$_3$H$_2$ increases rapidly because of the complete evaporation of c-C$_3$HD and c-C$_3$H$_2$ from the grain surface. Inside the inner 60 au, the gas-phase abundance of c-C$_3$HD decreases rapidly because of the slow deuteration process (i.e. unavailability of any forms of H$_2$D$^+$ to transfer deuteron). However, in this region, c-C$_3$H$_2$ can still survive with a very high abundance ($\sim 10^{-9}$) because of its efficient production via the reactions H + CH$_2$CCH → ortho-H$_2$ + c-C$_3$H$_2$ and H + CH$_2$CCH → para-H$_2$ + c-C$_3$H$_2$.

<table><caption>**Table 4.** Modelled and observed fractionations for c-C$_3$H$_2$.</caption><thead><tr><th>Radius</th><th>3000 au</th><th>4000 au</th><th>7000 au</th><th>Observed values<br>(at 3000 au)</th></tr></thead><tbody><tr><td>log<sub>10</sub> [c-C<sub>3</sub>HD]<sup>a</sup></td><td>-10.15</td><td>-10.07</td><td>-9.50</td><td>-11.3</td></tr><tr><td>log<sub>10</sub> [c-C<sub>3</sub>H<sub>2</sub>]<sup>a</sup></td><td>-8.37</td><td>-8.35</td><td>-8.33</td><td>-10.45<sup>c</sup></td></tr><tr><td>[c-C<sub>3</sub>HD/c-C<sub>3</sub>H<sub>2</sub>]<sup>b</sup></td><td>1.7%</td><td>2%</td><td>7%</td><td>14%</td></tr></tbody></table>

Notes. $^a$[X] = N(X)/N(H$_2$).

$^b$Deuterium fraction of c-C$_3$H$_2$.

$^c$Abundance of c-C$_3$H$_2$ by considering the same excitation temperature of c-C$_3$HD (see Section 2.3).

In Table 4, we give our modelled abundances and fractionation ratios at 3000, 4000 and 7000 au. The variation of predicted deuterium fractionation between 3000 and 4000 au is very small, whereas we observe a large variation at 7000 au (i.e. in the outer part of the envelope where the temperature is below 30 K). In the outer part of the envelope, H$_3^+$ reacts with HD, the major reservoir of D-atoms and D-atom is transferred from HD to H$_2$D$^+$ (ortho/para). As a result, H$_2$D$^+$ becomes very abundant in the outer part and serves as a primary species towards deuterium fractionation. Besides that, other abundant neutrals and important destruction partners of H$_3^+$ isotopologues, such as O and CO, deplete from the gas phase (e.g. because of the freeze out on to dust grains in cold and dense regions). As a consequence, in the outer part, the main formation reaction of c-C$_3$HD via c-C$_3$H$_2$D$^+$ + e$^-$ becomes very fast due to efficient deuteron transfer between para-H$_2$D$^+$ and c-C$_3$H$_2$ to form c-C$_3$H$_2$D$^+$, which results in an increase in the c-C$_3$HD abundance steeply beyond 5000 au. For c-C$_3$H$_2$, however, we have not seen similar behaviour because of its efficient destruction via various ion-molecular reactions with H$_3$O$^+$, para-H$_3^+$, HCO$^+$, ortho-H$_3^+$ and H$^+$. When we are going inside the envelope as compared to the outer part, deuteron transfer proceeds via other secondary species (e.g. at 2000 au, it is DCO$^+$ that originates from H$_2$D$^+$), which results in lower fractionation.

The deuterium fractionation predicted by our model at 3000 au (approximate size of the envelope) is 1.7 per cent (i.e. within a factor of 10 of what we obtain from the observations – about 14 per cent). The c-C$_3$HD modelled abundance is $7 \times 10^{-11}$ at 3000 au (see Table 4), which is within a factor of 10 of the observed value ($5 \times 10^{-12}$). For c-C$_3$H$_2$, our model overpredicts the observed value by two orders of magnitude. The discrepancy between the model and the observations could arise for a few possible reasons. First, in a highly centrally peaked source such as IRAS 16293, the physical structure is of utmost importance to determine the abundance profile of any species. Inhomogeneities in the protostellar envelope could affect the determination of the abundance profile. Although higher sensitivity observations and higher spectral resolution could help in this aspect (as the velocity profile will depend on the physical structure), the real solution would be higher spatial resolution (interferometric) observations. Secondly, it seems that current astrochemical models overpredict the abundance of c-C$_3$H$_2$ in different types of sources; see Agúndez & Wakelam (2013) for a review, and also Sipilä, Spezzano & Caselli (2016). Additional chemical studies of carbon chain productions would have to be carried out in order to reproduce the observed abundances.

## ACKNOWLEDGEMENTS

This work is based on observations carried out with the IRAM 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). LM, PG and VW are grateful for
ERC starting grant (3DICE, grant agreement 336474) that provided funding during this work. PG's postdoctoral position is funded by the INSU/CNRS. VW and PG acknowledges the French program 'Physique et Chimie du Milieu Interstellaire' (PCMI) funded by the Conseil National de la Recherche Scientifique (CNRS) and Centre National d'Etudes Spatiales (CNES). We would like to thank the anonymous referee for constructive comments that helped to improve the manuscript.

## REFERENCES

Agúndez M., Wakelam V., 2013, Chemical Reviews, 113, 8710

Aikawa Y., Wakelam V., Garrod R. T., Herbst E., 2008, ApJ, 674, 984

Bell M. B., Feldman P. A., Matthews H. E., Avery L. W., 1986, ApJ, 311, L89

Bell M. B., Avery L. W., Matthews H. E., Feldman P. A., Watson J. K. G., Madden S. C., Irvine W. M., 1988, ApJ, 326, 924

Bogey M., Demuynck C., Destombes J. L., Dubus H., 1987, Journal of Molecular Spectroscopy, 122, 313

Brown P. D., Millar T. J., 1989, MNRAS, 240, 25P

Caselli P., Walmsley C. M., Zucconi A., Tafalla M., Dore L., Myers P. C., 2002, ApJ, 565, 344

Caux E. et al., 2011, A&A, 532, A23

Ceccarelli C., Caselli P., Bockelée-Morvan D., Mousis O., Pizzarello S., Robert F., Semenov D., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. University of Arizona Press, Tucson, p. 859

Cox P., Guesten R., Henkel C., 1987, A&A, 181, L19

Crimier N., Ceccarelli C., Maret S., Bottinelli S., Caux E., Kahane C., Lis D. C., Olofsson J., 2010, A&A, 519, A65

Daniel F. et al., 2013, A&A, 560, A3

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Garrod R. T., Wakelam V., Herbst E., 2007, A&A, 467, 1103

Gerin M., Wootten H. A., Combes F., Boulanger F., Peters III W. L., Kuiper T. B. H., Encrenaz P. J., Bogey M., 1987, A&A, 173, L1

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5

Hasegawa T. I., Herbst E., 1993, MNRAS, 261, 83

Hincelin U., Wakelam V., Hersant F., Guilloteau S., Loison J. C., Honvault P., Troe J., 2011, A&A, 530, A61

Linsky J. L. et al., 2006, ApJ, 647, 1106

Lucas R., Liszt H. S., 2000, A&A, 358, 1069

Madden S. C., Irvine W. M., Swade D. A., Matthews H. E., Friberg P., 1989, AJ, 97, 1403

Majumdar L. et al., 2016a, MNRAS, 466, 4470

Majumdar L., Gratier P., Vidal T., Wakelam V., Loison J.-C., Hickson K. M., Caux E., 2016b, MNRAS, 458, 1859

Maret S., Hily-Blant P., Pety J., Bardeau S., Reynier E., 2011, A&A, 526, A47

Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350

Matthews H. E., Irvine W. M., 1985, ApJ, 298, L61

Müller H. S. P., Schlöder F., Stutzki J., Winnewisser G., 2005, Journal of Molecular Structure, 742, 215

Neufeld D. A., Wolfire M. G., Schilke P., 2005, ApJ, 628, 260

Ruaud M., Loison J. C., Hickson K. M., Gratier P., Hersant F., Wakelam V., 2015, MNRAS, 447, 4004

Sipilä O., Spezzano S., Caselli P., 2016, A&A, 591, L1

Spezzano S. et al., 2013, ApJ, 769, L19

Thaddeus P., Vrtilek J. M., Gottlieb C. A., 1985, ApJ, 299, L63

Vrtilek J. M., Gottlieb C. A., Thaddeus P., 1987, ApJ, 314, 716

Wakelam V., Vastel C., Aikawa Y., Coutens A., Bottinelli S., Caux E., 2014, MNRAS, 445, 2854

Wakelam V. et al., 2015, ApJS, 217, 20

Willacy K., Millar T. J., 1998, MNRAS, 298, 562

This paper has been typeset from a *TEX/LATEX* file prepared by the author.