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

TL;DR

This study reports the detection and analysis of singly deuterated cyclopropenylidene c-C3HD in the protostar IRAS 16293-2422, highlighting its potential as a tool for probing cold interstellar environments.

Contribution

First detection of c-C3HD in a protostellar environment and comparison with astrochemical models to understand deuterium fractionation.

Findings

c-C3HD can be used to study cold interstellar chemistry

Deuterium fraction in c-C3HD is approximately 14%

Results support the role of c-C3HD in astrochemical diagnostics

Abstract

Cyclic-C3HD (c-C3HD) is a singly deuterated isotopologue of c-C3H2, 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 setup using the EMIR heterodyne 3 mm receiver of the IRAM 30m telescope. We observed seven lines of c-C3HD and three lines of c-C3H2. Observed abundances are compared with astrochemical simulations using the NAUTILUS gas-grain chemical model. Our results clearly show that c-C3HD can be used as an important supplement for studying chemistry and physical conditions for cold environments. Assuming that the size of the protostellar envelope is 3000 AU and same excitation temperatures for both c-C3H2 and c-C3HD, we obtain a deuterium fraction of $14_{-3}^{+4}\%$.

Tables

Table 1. Table 1: Observed lines and spectroscopic parameters for c -C 3 HD and c -C 3 H 2

Species	Lines Ref.	Integrated flux	VLSR	FWHM	Observed	Aij	Eup	Quantum
		($K\mathrm{km}{s}^{-1}$)	($\mathrm{km}{s}^{-1}$)	($\mathrm{km}{s}^{-1}$)	Frequency	($s^{-1}$)	($K$)	numbers
					($\mathrm{MHz}$)
		$0.013\pm 0.009$	$4.3\pm 0.162$	$1.0\pm 0.70$	85643.318	$1.43\times {10}^{-5}$	26.6	$4_{3,2}\to 4_{2,3}$
		$0.028\pm 0.003$	$4.3\pm 0.064$	$1.3\pm 0.12$	102423.019	$1.53\times {10}^{-5}$	22.3	$4_{1,3}\to 4_{0,4}$
		$0.243\pm 0.024$	$4.2\pm 0.016$	$1.1\pm 0.02$	104187.126	$3.96\times {10}^{-5}$	10.8	$3_{0,3}\to 2_{1,2}$
c-C3HD	Bogey et al. (1987)	$0.047\pm 0.006$	$4.2\pm 0.031$	$0.8\pm 0.08$	104799.707	$7.29\times {10}^{-6}$	10.9	$3_{1,3}\to 2_{1,2}$
		$0.032\pm 0.005$	$4.5\pm 0.086$	$1.4\pm 0.19$	106256.108	$1.69\times {10}^{-5}$	22.5	$4_{2,3}\to 4_{1,4}$
		$0.054\pm 0.006$	$4.1\pm 0.025$	$0.9\pm 0.07$	106811.09	$7.87\times {10}^{-6}$	10.8	$3_{0,3}\to 2_{0,2}$
		$0.280\pm 0.028$	$4.1\pm 0.005$	$0.8\pm 0.01$	107423.671	$4.47\times {10}^{-5}$	10.9	$3_{1,3}\to 2_{0,2}$
		$0.125\pm 0.013$	$4.2\pm 0.018$	$1.4\pm 0.06$	84727.688	$1.04\times {10}^{-5}$	16.14	$3_{2,2}\to 3_{1,3}$
c-C3H2	Thaddeus et al. (1985)	$1.372\pm 0.137$	$4.2\pm 0.010$	$1.4\pm 0.006$	85338.894	$2.32\times {10}^{-5}$	6.45	$2_{1,2}\to 1_{0,1}$
		$0.173\pm 0.017$	$4.3\pm 0.028$	$1.9\pm 0.05$	85656.431	$1.52\times {10}^{-5}$	29.07	$4_{3,2}\to 4_{2,3}$

Table 2. Table 2: Priors distribution functions for the parameters used in the bayesian approach.

Parameter	Distribution
${\log }_{10}N$ (${\mathrm{cm}}^{-2}$)	$\mathrm{Uniform}(8,22)$
${\log }_{10}{T}_{\mathrm{ex}}$ ($K$)	$\mathrm{Uniform}({\log }_{10}3,{\log }_{10}200)$
$V$ ($\mathrm{km}{s}^{-1}$)	$\mathrm{Uniform}(2.8,4.8)$
$\Delta V$ ($\mathrm{km}{s}^{-1}$)	$\mathrm{Uniform}(0.25,10)$
${\log }_{10}{\sigma }_{\mathrm{add}}$ ($K$)	$\mathrm{Uniform}(-3,1)$

Table 3. Table 3: Point estimates of the posterior distribution function corresponding to the median and one sigma uncertainty.

Parameter	c-C3HD	c-C3H2
${\log }_{10}N$ (${\mathrm{cm}}^{-2}$)	$12.83\pm 0.05$	$13.42\pm 0.01$
${\log }_{10}{T}_{\mathrm{ex}}$ ($K$)	$0.80\pm 0.05$	$0.95\pm 0.03$
${\log }_{10}\Delta V$ ($\mathrm{km}{s}^{-1}$)	-$0.09\pm 0.02$	$0.11\pm 0.01$
${\log }_{10}{\left[X\right]}^a$	$-11.3\pm 0.05$	$-10.7\pm 0.01$
Notes: a [X] = N(X)/N($H_2$) and b value with $95\%$ confidence.
Here quoted uncertainties are statistical only, without the
10% calibration error.

Table 4. Table 4: Modeled and observed fractionations for c -C 3 H 2 .

Radius	3000 AU	4000 AU	7000 AU	Observed values
				( at 3000 AU)
${\log }_{10}{[c-C_3\mathrm{HD}]}^a$	-10.15	-10.07	-9.50	-11.3
${\log }_{10}{[c-C_3{H}_2]}^a$	-8.37	-8.35	-8.33	$\text{-}{10.45}^c$
${\left[c-C_3\mathrm{HD}/c-C_3{H}_2\right]}^b$	1.7%	2%	7%	14%
Notes: a [X] = N(X)/N($H_2$)
Notes: b Deuterium fraction of c-C3H2
Notes: c Abundance of c-C3H2 by considering the same excitation temperature of c-C3HD (see Section 2.3)

Equations

$$\begin{array}{r}{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 }],\end{array}$$

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

$$\begin{array}{rl}\frac{\mathrm{err}{\mathrm{N}}_2}{{\mathrm{N}}_2}& =\sqrt{(\frac{\mathrm{errstat}{\mathrm{N}}_2}{{\mathrm{N}}_2}{)}^2+(\mathrm{errcal}{)}^2}\\ & =\sqrt{0.02^2+0.1^2}=0.10\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{err}{\mathrm{N}}_1}{{\mathrm{N}}_1}& =\sqrt{(\frac{\mathrm{errstat}{\mathrm{N}}_1}{{\mathrm{N}}_1}{)}^2+(\mathrm{errcal}{)}^2}\\ & =\sqrt{0.1^2+0.1^2}=0.14\end{array}$$

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

$$\begin{array}{rl}{\mathrm{R}}_{\min }& ={\mathrm{N}(c\text{mbox}-{\mathrm{C}}_3\mathrm{HD})}_{\min }/{\mathrm{N}(c\text{mbox}-{\mathrm{C}}_3{\mathrm{H}}_2)}_{\max }\\ & =6.7\times (1-0.14)\times 10^{12}/4.7\times (1+0.1)\times 10^{13}\\ & =11\%\end{array}$$

Full text

††pubyear: 2015

1 Introduction

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

c-C3H2 shows strong rotational transitions in millimetre wave bands due to 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 the discussion). c-C3H2 is then often used as a template to study deuteration along different sources since the rotational transitions of its isotopic species are also very intense. c-C3H2 can exist in singly (c-C3HD) and doubly (c-C3D2) deuterated forms. c-C3HD 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%) of the $2_{1,2}\rightarrow 1_{0,1}$ lines of c-C3HD to c-C3H2 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 (TROT= 4-7 K) and thus c-C3HD samples the coldest part of the interstellar medium (ISM). c-C3D2 was also detected recently by Spezzano et al. (2013) toward the starless cores TMC-1C and L1544. They measured the abundance of c-C3D2 with respect to the c-C3H2 and found it to be 0.4%-0.8% in TMC-1C and 1.2%-2.1% in L1544. This clearly shows that the abundance of c-C3D2 is much lower than that of c-C3HD and thus supports the idea that singly deuterated species forms more easily in the ISM as compared to the multiply deuterated species (Brown & Millar (1989), Willacy & Millar (1998)). This makes c-C3HD a better candidate to study the chemistry and the physical conditions of cold environments compared to c-C3D2. The main goal of the present paper is to report the detection of c-C3HD in the low mass protostar IRAS 16293-2422 (hereafter IRAS 16293) and study deuterium fractionation of c-C3H2 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 and finally results are discussed in the last Section.

2 Observations and data reduction

2.1 Observations

The observations were obtained at the IRAM 30m telescope during the period August 18-23, 2015. Overall, the weather condition was like an average summer (a median value of 4-6 mm water vapour). We performed our observations by using the EMIR heterodyne 3 mm receiver tuned at a frequency of 89.98 GHz in the Lower Inner sideband. This receiver was followed by a Fourier Transform Spectrometer in its 195$\mathrm{\,kHz}$ resolution mode. The observed spectrum was composed of two regions: one from 84.4 GHz to 92.3 GHz and another one from 101.6 GHz to 107.9 GHz.

The observations were made towards the midway point between sources A and B of IRAS 16293 at $\mathrm{\alpha_{2000}=16^{h}32^{m}22.75^{s}}$, $\mathrm{\delta_{2000}=-24^{\circ}28^{{}^{\prime}}34.2^{{}^{\prime\prime}}}$. The A and B components, separated by $5.5^{\prime\prime}$, 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 seconds and a throw of $90^{\prime\prime}$ 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^{\prime\prime}$).

2.2 Results

2.2.1 c-C3HD and c-C3H2 line properties

The data were reduced and analysed using the CLASS software from the GILDAS111https://www.iram.fr/IRAMFR/GILDAS/ package. We made Gaussian fits to the detected lines following a local low (0 or 1) order polynomial baseline subtraction. Table 1 shows the result of these fits for the 7 observed lines of c-C3HD and 3 observed lines of c-C3H2. In Table 1, uncertainties on integrated flux are the Gaussian fit uncertainties (with an added 10% calibration error), uncertainties on VLSR are the quadratic sum of the gaussian fit statistical uncertainties and the frequency uncertainties from the spectroscopic catalog, uncertainties on FWHM are the gaussian fit uncertainties only. All the observed lines are single component features with the mean LSR velocity of $\sim$ 4.2 km/s and the mean FWHM of $\sim$ 1 km/s. Caux et al. (2011) identified four types of kinematical behaviours for different species based on the different FWHM and VLSR distributions. We found that c-C3HD and c-C3H2 belong to the type I (i.e. FWHM$\leq$2.5 km/s, VLSR$\sim$ 4 km/s, Eup$\sim$ 0-50 K), which corresponds to species abundant in the cold envelope of IRAS 16293.

2.2.2 Radiative transfer modelling for c-C3HD and c-C3H2

Majumdar et al. (2016b) used a bayesian model to recover the distribution of parameters which best agree with the observed line intensities for $\mathrm{CH_{3}SH}$ in the same source. Here, we use the same model with the exception that not only the integrated intensities but also the full modelled spectra are compared to the observed ones.

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-C3HD and c-C3H2. The input parameters in this radiative transfer model are the species column density, the line width, the excitation temperature, the source size along with an accurate spectroscopic catalog. The spectroscopic catalogs for c-C3HD and c-C3H2 were retrieved from the 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^{\prime\prime}$ ($\sim$ 3000 AU, typical size of the protostellar envelope by Caux et al. (2011)) since 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 EMCEE222https://github.com/dfm/emcee (Foreman-Mackey et al., 2013) of the Affine Invariant Ensemble Monte Carlo Markov Chain approach (Goodman & Weare, 2010). Sixty walkers are initialised in a small 6 dimensional sphere in the center of the parameter space. The chains are evolved for a burn-in sequence of 3000 steps after which the convergence is checked by examining a plot of the running mean of each parameter.

Figure 1 shows the comparison of observed spectra and modelled spectra. Figures 2 and 3 show the 1D and 2D histograms of the posterior probability distribution function for c-C3HD and c-C3H2 respectively. In both cases, all parameters are well defined. Table 3 summarises the one point statistics for the marginalised posterior distributions of parameters. The uncertainties reported in Table 3 are 1$\sigma$ symmetric error bars. There is about a factor 2 difference in line width between c-C3HD and c-C3H2 (shown in Table 3), this cannot arise from opacity effects as all the lines are optically thin except c-C3H2 line at 85.338 GHz. One possibility is that c-C3H2 comes both from the central hotter region and the outer envelope while c-C3HD only comes from the outer cold envelope. This is further strengthen by the chemical modelling which shows such a pattern in Figure 4. Thus, 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-C3H2

The observed main beam temperature $T_{mb}$, for a beam of size $\theta_{t}$ of a gaussian source of size $\theta_{s}$ is

[TABLE]

where, $\eta$ is the beam dilution factor which is identical for two molecules since 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 the optically thin conditions, we can assume that $\tau\propto N$. Under this assumption, $T_{mb}\propto[J(T_{ex})-J(T_{CMB})]N$.

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 doing the independent analysis. We need more lines, 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 modeling of isotopologues (See for instance Caselli et al. (2002); Daniel et al. (2013)).

Now we can represent the column density (N2) of c-C3H2 at the excitation temperature (6.3 K) of c-C3HD as:

[TABLE]

Where, N1 is the column density of c-C3H2 at 8.9 K. At frequency 85 GHz, $J(T_{ex1})=7.02$ for $T_{ex1}$ of 8.9 K; $J(T_{ex2})=4.46$ for $T_{ex2}$ of 6.3 K [$J(T_{CMB})=1.17$ at $T_{CMB}=2.73$ K], we derive a column density of $\rm{N({\it c}}$-$\rm{C_{3}H_{2})}$=$N_{2}=1.78N_{1}=4.7\times 10^{13}\mathrm{\,cm^{-2}}$. Thus, we obtain a deuteration fraction $\rm{N({\it c}}$-$\rm{C_{3}HD)}$/$\rm{N({\it c}}$-$\rm{C_{3}H_{2})}$=$6.7\times 10^{12}/4.7\times 10^{13}=14\%$.

If we define the error on column densities of c-C3H2 and c-C3HD by errN2 and errN$1$, statistical error on column densities of c-C3H2 and c-C3HD by errstatN2 and errstatN1 and relative calibration error by errcal (10% in both cases), then we can write:

[TABLE]

[TABLE]

For N(c-C3HD)/N(c-C3H2)=R, then

[TABLE]

[TABLE]

Thus, error on the deuterium fraction can be written as N(c-C3HD)/N(c-C3H2)=$14_{-3}^{+4}\%$. Bell et al. (1988) also measured deuterium fractionation of c-C3H2 in the order of 5 to 15% in TMC-1.

3 1D modelling of the protostellar envelope

3.1 The NAUTILUS chemical model with deuteration

To model the chemistry of c-C3HD and c-C3H2 in IRAS 16293, we used the same approach as Majumdar et al. (2016b). For our simulation, we have used the 2 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. submitted). 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 2 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 onto grain surfaces, diffusion of the accreted species, reaction at the grain surface, and finally by 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 et al. (2007). We adopt the similar initial elemental abundances 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 et al., 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 H2 and HD forms respectively. For our standard model, we have used a C/O ratio of 0.7 and an ortho-to-para H2 ratio of 3.

In our current model, all the variables relative to H2 has been modified in terms of nuclear spin states ortho-H2 and para-H2. Since, it has been known from long time that H2, D2, H3+, H2D+, D2H+ and D3+ along with their spin isomers are the main species that dictate deuterium fractionation at low temperature (see Ceccarelli et al., 2014). Detailed description and benchmarking of our deuterated network with spin chemistry is presented in Majumdar et al. (2016a). This network is available on the KIDA333http://kida.obs.u-bordeaux1.fr/ website. 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 1D physical structure

To follow the deuterium fractionation of c-C3H2 in IRAS 16293, we have used the same 1D physical structure as in Aikawa et al. (2008); Wakelam et al. (2014); Majumdar et al. (2016b). This physical structure was based on the nongray radiation hydrodynamic model by Masunaga & Inutsuka (2000) to follow the core evolution from pre-stellar core to protostellar core. It starts from a hydrostatic prestellar core with central density $n$(H2) $\sim 3\times 10^{4}$ cm*-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. Initial temperature for the core is around 7 K at the center 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, core 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 rising in temperature in the central region. Contraction then decelerated due to increase in the gas pressure which eventually makes the first hydrostatic core, known the ‘first core’ at the center. When the core center reaches very high density ($10^{7}$ cm*-3*) and high temperature (2000 K), the hydrostatic core becomes unstable due to H2 dissociation and starts to collapse again. This collapse is referred as the ‘second collapse’. Within a short period of time, the dissociation degree approaches unity at the center due to 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 ‘protostellar core’. In the model, the initial prestellar core evolves to the protostellar core in $2.5\times 10^{5}$ yr. When 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 discussions

Here, we present the predicted abundances of c-C3HD and c-C3H2 in the protostellar envelope using the physical and chemical models previously described, and compare with our observations. Figure 4 shows the computed abundances of c-C3HD and *c-*C3H2, 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 since 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 multi-wavelength dust and molecular observations (see Wakelam et al. (2014) for the discussion).

Abundance profiles of c-C3HD and c-C3H2 predicted by our model can be divided into three regions. 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-C3HD and c-C3H2 decrease toward 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-C3HD forms mainly by the dissociative recombination of c-C3H2D+. c-C3H2D+ is produced mainly from the deuteron transfer from para-H2D+ and DCO+ to c-C3H2. c-C3H2D+ is also partly produced by c-C3H2 + H2DO+ and para-H2 + C3D+ reactions. So in the cold outer envelope, deuteration of c-C3HD is mainly controlled by para-H2D+ and DCO+. In these region, c-C3H2 also mainly forms from the dissociative recombination reaction of c-C3H3+. From 2000 AU to 200 AU (where temperature started to increase from 30 K to 50 K), c-C3HD mainly forms from the neutral-neutral reaction CH + C2HD $\rightarrow$ H + c-C3HD whereas c-C3H2 mainly forms from the H + CH2CCH $\rightarrow$ para-H2 + c-C3H2 reaction.

Between 200 and 100 AU, the gas phase abundance of both c-C3HD and c-C3H2 increases rapidly due to the complete evaporation of c-C3HD and c-C3H2 from the grain surface. Inside the inner 60 AU, the gas phase abundance of c-C3HD decreases rapidly due to slow deuteration process, i.e. unavailability of any forms of H2D+ to transfer deuteron. But in this region, c-C3H2 can still survive with a very high abundance ($\sim 10^{-9}$) due to its efficient production via the reactions H + CH2CCH $\rightarrow$ ortho-H2 + c-C3H2 and H + CH2CCH $\rightarrow$ para-H2 + c-C3H2.

In Table 4, we give our modelled abundances and fractionation ratios at 3000 AU, 4000 AU and 7000 AU. The variation of predicted deuterium fractionation between 3000 AU 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, H3+ reacts with HD, the major reservoir of D-atoms and D-atom is transferred from HD to H2D+ (ortho/para). As a result, H2D+ 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 H3+ isotopologues, such as O and CO, deplete from the gas phase (for example because of the freeze out onto dust grains in cold and dense regions). As a consequence, in the outer part, the main formation reaction of c-C3HD via c-C3H2D+ + e- becomes very fast due to efficient deuteron transfer between para-H2D+ and c-C3H2 to form c-C3H2D+ which results into an increase in the c-C3HD abundance steeply beyond 5000 AU. For c-C3H2, however we have not seen similar behavior due to its efficient destruction via various ion-molecular reactions with H3O+, para-H3+, HCO+, ortho-H3+, and H+. When we are going inside the envelope as compared to outer part, deuteron transfer proceeds via other secondary species (for example at 2000 AU, it is DCO+ which originates from H2D+) which results into lower fractionation.

The deuterium fractionation predicted by our model at 3000 AU (approximate size of the envelope), is of 1.7%, i.e. within a factor of 10 than what we obtain from the observations (about 14%). The c-C3HD 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-C3H2, our model over predicts the observed one by two orders of magnitude. The discrepancy between the model and the observations could arise from few possible reasons. First, in a highly centrally peaked source like IRAS 16293, the physical structure is of utmost importance to determine the abundance profile of any species. Inhomogeneities in the protostellar envelop could affect the determination of the abundance profile. Although, higher sensitivity observations and higher spectral resolution could help on this aspect (since the velocity profile will depend on the physical structure), the real solution would be higher spatial resolution (interferometric) observations. Second, it seems that current astrochemical models over predicts the abundance of c-C3H2 in different types of sources (see Agúndez & Wakelam (2013) for the review; Sipilä et al. (2016)). Additional chemical studies of the carbon chain productions would have to be done to reproduce the observed abundances.

Acknowledgements

Based on observations carried out with the IRAM 30m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). LM, PG, VW thanks ERC starting grant (3DICE, grant agreement 336474) for funding during this work. PG postdoctoral position is funded by the INSU/CNRS. VW acknowledge the CNRS programme PCMI for funding of their research. We would like to thank the anonymous referee for constructive comments that helped to improve the manuscript.