H₂CO ortho-to-para ratio in the protoplanetary disk HD 163296

We observed the $\mathrm{H_{2}CO}$ $4_{04}-3_{03}$ line at 290.623 GHz with ALMA in May 2014 during Cycle 1 (project ADS/JAO.ALMA#2012.1.00681; PI: C. Qi.) with 32 antennas. The Band 7 observations have baselines lengths spanning between 27 and 650 m, and a total on-source time of 45.5 min. The correlator was configured with four spectral windows, covering the N₂H⁺ $J=3-2$, DCO⁺ $J=4-3$ and the $\mathrm{H_{2}CO}$ $4_{04}-3_{03}$ lines, as well as a window dedicated to continuum. A detailed description of the observations and the results of the N₂H⁺ $J=3-2$ and DCO⁺ $J=4-3$ lines are given in Qi et al. (2015). The $\mathrm{H_{2}CO}$ line was covered with a spectral window of 58 MHz bandwidth and 15 kHz channel width. To calibrate the amplitude and phase temporal variations, the quasar J1733-1304 was observed.

We used the standard calibration performed by the ALMA staff, and further self-calibrated the data in phase and amplitude to improve the signal-to-noise ratio. The continuum subtracted visibilities were cleaned in CASA 4.7, using the CLEAN algorithm with Briggs weighting of 2. The data were regridded to a spectral resolution of 0.6$\mathrm{\,km\,s^{-1}}$. We applied a Gaussian taper of $1^{\prime\prime}$ to improve the signal-to-noise. In addition, we created a mask to help the cleaning process, by selecting regions in each channel with line emission that is consistent with the expected Keplerian rotation of the disk. The resulting noise properties and beam information are summarized in Table 1.

We observed five $\mathrm{H_{2}CO}$ lines with the SMA. The $\mathrm{p-H_{2}CO}$ lines at 218 GHz and 290 GHz, and the $\mathrm{o-H_{2}CO}$ line at 300 GHz were observed in March 2014 with 7 antennas in the sub-compact (SUB) array configuration. The 218 GHz line was observed in one spectral setting, with a spectral window of 104 MHz and 406 kHz channel width. The 290 and 300 GHz lines were covered in another setting, with two spectral windows of 104 MHz and 406 kHz channel width each. The quasars NRAO 530 and $1744-312$ were observed to calibrate amplitude and phase variations, and the frequency bandpass was calibrated using observations of 3C279. Titan and Mars were used as flux calibrators. We also include observations of the $\mathrm{o-H_{2}CO}$ $4_{14}-3_{13}$ and $3_{12}-2_{11}$ lines at 281.527 and 225.698 GHz, respectively, which have been presented in Qi et al. (2013). A detailed description of the observations can be found in Qi et al. (2013). For consistency with the other $\mathrm{H_{2}CO}$ lines, we present an independent imaging in this paper.

The SMA data were first calibrated with the MIR software package using standard procedures¹¹1https://www.cfa.harvard.edu/~cqi/mircook.html. The calibrated data were then exported to CASA and cleaned in the same way as the ALMA line. The data were regridded to a common spectral resolution of 0.6$\mathrm{\,km\,s^{-1}}$. A robust parameter of 2 was used to optimize the sensitivity. To help the cleaning process, we created circular masks for each line, that was the same for each channel. The radius of the mask was chosen to cover the line emission in all channels, and varied for the different $\mathrm{H_{2}CO}$ lines. The resulting rms noise levels and observational parameters are summarized in Table 1.

Fig. 1 shows the ALMA observations of the $\mathrm{p-H_{2}CO}$ $4_{03}-3_{03}$ line. The two left panels show the velocity integrated map, integrated over the full velocity range (left) and over the blue- and red-shifted part of the line (right). The third panel in Fig. 1 shows the azimuthally-averaged deprojected profile of the $\mathrm{p-H_{2}CO}$ $4_{03}-3_{03}$ line emission. It was obtained by extracting elliptical anulii in the unclipped moment-zero map, taking into account the inclination and position angle of the disk. The right-most panel shows the disk-integrated spectra, which was obtained using an elliptical mask. The integrated spectra shows a double-peaked profile, typical of Keplerian rotation of an inclined disk. The emission presents a ring-like spatial distribution, with a depletion of $\mathrm{H_{2}CO}$ emission at the disk center, a peak of emission at intermediate radii ($\sim 100$ au), and either a second but weaker peak or a plateau in the outer disk near $\sim 400$ au. This second bump is located farther out in the disk than the milimeter dust edge near $\sim 250$ au (Isella et al., 2016). Carney et al. (2017) presented higher-angular resolution observations (0.5”) of the lower energy ($E_{u}=21$ K) $3_{03}-2_{02}$ line, which has a similar radial profile.

The observed $\mathrm{H_{2}CO}$ distribution is not unique to the HD 163296 disk. Spatially resolved observations of the low-energy ($E_{u}\sim 30$ K) $\mathrm{H_{2}CO}$ $3-2$ line in TW Hya also show a depression of the emission at the disk center and a bump in the outer disk (Öberg et al., 2017). In contrast, observations of the higher-energy ($E_{u}\sim 60$ K) $5-4$ line in DM Tau and TW Hya present centrally peaked emission profiles (Loomis et al., 2015; Öberg et al., 2017). Such different morphologies of the $3-2$ and $5-4$ lines are well explained by the presence of hot $\mathrm{H_{2}CO}$ close to the star followed by a lack of $\mathrm{H_{2}CO}$ until cold $\mathrm{H_{2}CO}$ formation becomes possible in the outer disk. In summary, the observed $\mathrm{H_{2}CO}$ morphology toward HD 163296 is consistent with previous $\mathrm{H_{2}CO}$ disk observations, suggesting that it is a good object for a case study on $\mathrm{H_{2}CO}$ chemistry in disks.

The SMA observations are shown in Fig. 2. The top row shows the velocity integrated maps of the five $\mathrm{H_{2}CO}$ lines, while the red- and blue-shifted parts of the emission are shown in the second row. The bottom row shows the disk-integrated spectra, obtained using a circular mask with a radius of 6”. Although the emission looks washed-out in the moment-zero maps for some of the lines, visual inspection of the individual channels show line emission consistent with the expected Keplerian rotation of the disk. The two higher energy $\mathrm{H_{2}CO}$ lines ($E_{u}>45$ K) are centrally peaked. The peak emission of the remaining lines with $E_{u}<40$ K are not centered, suggesting ring-like emission patterns or a drop of the emission in the inner disk consistent with the higher-angular resolution ALMA observations.

We extracted disk-integrated fluxes using a circular mask with a radius of 6”. The fluxes are reported in Table 2. The uncertainty in the flux was estimated by taking the standard deviation of 500 simulated flux measurements from regions free of signal, using the same circular mask but centered at random positions in the map. The fluxes of the 290 GHz line observed with ALMA and SMA are consistent within the uncertainties. We use the disk-integrated fluxes to estimate an excitation temperature for $\mathrm{H_{2}CO}$.

The critical densities of the observed $\mathrm{H_{2}CO}$ lines, defined as the density at which the spontaneous radiative de-excitation rate ($A_{ul}$) is equal to the collisional de-excitation rate ($\gamma_{ul}$), are given in Table 2 for a temperature of 20 K. For a higher temperature of 60 K, the critical densities increase by less than 20%. When the gas density is comparable or larger than the critical density, the emission is thermalized and the excitation temperature is equal to the gas kinetic temperature. When the gas density is lower, then the emission becomes sub-thermally excited and $T_{\mathrm{ex}}<T_{\mathrm{kin}}$. Given the high ($>10^{6}~\rm{cm}^{-3}$) gas densities in protoplanetary disks, in particular closer to the midplane (see Fig. 4), we expect the bulk of the emission to be thermalized.

Using the three ortho lines we constructed a rotational diagram (Fig. 3). We obtain a rotational temperature of $24\pm 14$ K. This temperature is consistent with the temperature of 24 K derived by Carney et al. (2017), from the $\mathrm{p-H_{2}CO}$ $3_{03}-2_{02}$, $3_{22}-2_{21}$ and $3_{21}-2_{20}$ lines, using a match-filter technique to extract the fluxes of the two weaker lines. The derived rotational temperature also matches the CO freeze-out temperature found by Qi et al. (2015) for this disk, using N₂H⁺ and C¹⁸O observations. This suggest that the $\mathrm{H_{2}CO}$ emission at the large scales, to which the SMA is sensitive, arises predominantly from the cold layers close to the midplane, where CO starts to freeze-out onto dust grains and can be hydrogenated to form $\mathrm{H_{2}CO}$. In the next section, we use this measured rotational temperature to inform our model on the distribution of $\mathrm{H_{2}CO}$ in the disk.

Based on chemical model predictions of disks, the distribution of $\mathrm{H_{2}CO}$ can be separated in two main components: 1) A warm component towards the center of the disk, that is mainly produced by gas-phase chemistry (e.g., Loomis et al., 2015), and 2) a colder component located in the layers above the midplane in the outer disk, where the $\mathrm{H_{2}CO}$ formation is dominated by grain surface chemistry (e.g., Öberg et al., 2017). In the midplane, $\mathrm{H_{2}CO}$ remains frozen on the surfaces of dust grains. Following this, we parameterize the abundance structure of $\mathrm{H_{2}CO}$ assuming the two components displayed in Fig. 4. The first component is located in the inner disk, between $R_{1}$ and $R_{2}$, with constant abundance $X_{1}$. The second component is located between $R_{2}$ and $R_{out}$, and has an abundance described by a power-law, $X=X_{2}(r/100~\mathrm{au})^{\alpha}$. Both components have lower and upper boundaries, $(z/r)_{l}$ and $(z/r)_{u}$.

Given the uncertainties when modeling only one transition, we fixed some of the parameters. The fixed parameters are $R_{1},R_{2},R_{out},(z/r)_{l}$ and $(z/r)_{u}$; and only $X_{1},X_{2},$ and $\alpha$ are left free in the line fitting. The adopted values in the model are given in the left column of Table 3. The values for $R_{1}$ and $R_{2}$ are chosen based on predictions from chemical models of this disk (Cleeves et al., private communication). $R_{out}$ is fixed to 600 au, corresponding to the extent of the $\mathrm{H_{2}CO}$ emission. We note that due to the low-angular resolution of the SMA observations, the inferred disk-averaged o/p ratio (see section 4.2) is not very sensitive to variations in the position of the inner component (i.e., $R_{1}$ and $R_{2}$). We consider three different values (0.15, 0.20 and 0.25) for the lower boundary, $(z/r)_{l}$, to investigate the impact of the adopted height of the $\mathrm{H_{2}CO}$ disk emission layer on the derived o/p ratio. The upper boundary, $(z/r)_{u}$, is fixed to 0.5; we tried different values, and found that the results were not sensitive to this paramter.

We use a Bayesian approach to find the best fit of the model to the ALMA visibilities. For this, we first create a synthetic observation of the $\mathrm{p-H_{2}CO}$ line emission using the RADMC-3D package (Dullemond et al., 2012) to compute the level populations, assuming the gas is under local thermal equilibrium (LTE). We then compute the Fourier Transform of the modeled line emission, using the vis-sample Python package²²2https://pypi.python.org/pypi/vissample that re-projects the modeled visibilities on the observed $u-v$ points. We then compute the weighted difference between model and observations, for the real and imaginary parts of the complex visibility, which gives the likelihood of the model. This process is repeated by sampling the posterior distribution with the MCMC method implemented in the emcee package by Foreman-Mackey et al. (2013).

The best fit values of the model parameters are listed in the right columns in Table 3, for the three different lower boundaries we consider. Because we only model one $\mathrm{p-H_{2}CO}$ line, the abundance we derive corresponds to that of $\mathrm{p-H_{2}CO}$. The abundances listed in Table 3, however, correspond to the total $\mathrm{H_{2}CO}$ abundance, including both ortho and para, and are computed assuming a uniform o/p ratio of 3. We allow this parameter to vary in the next section. The best-fit is found for the model with $(z/r)_{l}=0.15$. Fig. 5 shows channel maps of the $\mathrm{H_{2}CO}$ line observed with ALMA (top row), together with the best fit model (middle row), and the residuals (bottom row). The inferred abundances in the inner $\sim 200$ au change considerably (by more than one order of magnitude) depending on the adopted value for the lower boundary. The $\mathrm{H_{2}CO}$ abundance in the outer disk, given by the second component, change by a factor of $\sim 3-7$ (at $300-500$ au) within the three sets of models.

Motivated by the bump observed in the radial profile of the $\mathrm{H_{2}CO}$ ALMA line (see Fig. 1), we also explored the presence of a third component in the outer disk that is closer to the midplane. For this, we ran another set of models including a third component of constant abundance between 350 and 600 au, and upper boundary equal to the lower boundary of the second component (i.e., $0<z/r<(z/r)_{l}$). The best-fit abundance of this third component remained low compared to the two other components no matter the adopted $(z/r)_{l}$ value, suggesting that the contribution of such a third cold component to the overall $\mathrm{H_{2}CO}$ emission is small. We note that these models improved slightly the fit of the ALMA $\mathrm{H_{2}CO}$ line but resulted in similar results for the SMA observations.

We model the line emission of the two para and three ortho lines of $\mathrm{H_{2}CO}$ in the same way we did for the ALMA observations. Given the lower angular resolution of the SMA observations, we keep the $\mathrm{H_{2}CO}$ abundance structure fixed, and vary only the o/p ratio. We assume a uniform o/p ratio across the disk. In order to take into account for differences in the absolute fluxes between the ALMA and SMA data due to absolute flux calibration uncertainites of 10-20%, we introduce a scale factor for the total $\mathrm{H_{2}CO}$ abundance in the disk. That is, the total $\mathrm{H_{2}CO}$ abundance derived in the previous section is multiplied by this scale factor, which is a free parameter in the SMA line fitting.

The resulting best-fit parameters are shown in the bottom rows of Table 3. Fig. 6 shows the posterior probability distribution for the $\mathrm{H_{2}CO}$ disk-averaged o/p ratio. The scale factor is $1.32-1.76$ (depending on $(z/r)_{l}$), and increases somewhat with increasing lower boundary. These values are of the expected magnitude when considering the differences in measured fluxes with ALMA and SMA (see section 3.2).

The inferred o/p ratio changes depending on the assumed lower boundary but remains consistently lower than 3, the high-temperature value. The highest o/p value (2.5) is obtained when $\mathrm{H_{2}CO}$ is located closer to the midplane, and lower values (2.2 and 2.0) are obtained when $\mathrm{H_{2}CO}$ is closer to the disk surface. The possible values for the $\mathrm{H_{2}CO}$ o/p ratio are $1.8-2.8$ (within 1$\sigma$), corresponding to spin temperatures of $11-22$ K. Fig. 3 shows the rotational diagram of $\mathrm{o-H_{2}CO}$, and compares the observed and best-fit model integrated fluxes for the three $(z/r)_{l}$ choices. The best-fit is obtained for $(z/r)_{l}=0.15$, which gives an o/p ratio of $2.50^{+0.31}_{-0.27}$, corresponding to $T_{\mathrm{spin}}=17^{+6}_{-3}$ K. This is the only model resulting in consistent spin and rotational temperatures. However, we cannot rule out any of the vertical distribution choices based on the SMA data.

We note that this corresponds to a disk-averaged o/p ratio. Most likely, this value changes across the disk, both radially and vertically. Our inferred o/p value thus suggest that a considerable amount of the gas in the HD 63296 disk has an $\mathrm{H_{2}CO}$ o/p ratio that is lower than 3. Higher-angular resolution observations are needed to disentangle regions in the disk with low o/p ratios from regions presenting a high-temperature statistical value.

The $\mathrm{H_{2}CO}$ o/p ratio has been measured in a variety of interstellar environments and one comet. Fig. 7 presents a compilation from the literature of measured $\mathrm{H_{2}CO}$ o/p ratios. In the subsequent panels we present the same data on $\mathrm{H_{2}O}$ and $\mathrm{NH_{3}}$. We highlight with red edges those sources where the o/p ratio has been measured for $\mathrm{H_{2}CO}$ and another species. The solid curves in Fig. 7 correspond to the predicted o/p curves given by the Boltzmann distribution of the ortho and para energy levels. The observations are also listed in Table 4.

A range of $\mathrm{H_{2}CO}$ o/p ratios is found in the ISM. Low-temperature ($T_{\mathrm{spin}}<30$ K) values have been observed in different environments, such as the envelopes around low-mass stars (Jørgensen et al., 2005), outflows (Dickens & Irvine, 1999), and the Horsehead PDR (Guzmán et al., 2011). $\mathrm{H_{2}CO}$ o/p ratio consistent with the high-temperature value of 3 have been observed towards quiescent cores, such as the cold ($\sim 20$ K) UV-shielded core located just behind the Horsehead PDR (Guzmán et al., 2011), and the Orion-Bar PDR (Cuadrado et al., 2017). Because the Orion-Bar ($T_{\mathrm{dust}}\sim 150$ K) is a much stronger UV-illuminated PDR, dust grains at the cloud edge are free from ice mantles, in contrast to the Horsehead PDR where icy mantles can exist ($T_{\mathrm{dust}}\sim 20-30$ K). Cuadrado et al. (2017) have thus suggested that $\mathrm{H_{2}CO}$ forms predominantly in the gas-phase in the Orion-Bar, while it forms on ices and is later non-thermally desorbed into the gas-phase in the Horsehad PDR. However, recent results from laboratory experiments cast doubt on this interpretation (see section 5.2 for a discussion).

The $\mathrm{H_{2}CO}$ o/p ratio has been measured in one comet, 103P/Hartley 2, where a value of $2.1\pm 0.6$ was found, corresponding to $T_{\mathrm{spin}}\sim 13.5^{+6.7}_{-3.3}$ K (Gicquel et al., 2014). This o/p value is consistent with our disk measurement of $1.8-2.8$, but the spin temperature is lower compared to what has been deduced from a sample of $\mathrm{NH_{3}}$ and $\mathrm{H_{2}O}$ o/p ratio comet measurements; Shinnaka et al. (2016) measured the $\mathrm{NH_{3}}$ o/p ratio in 26 comets, and found values of $1.1-1.2$, corresponding to a spin temperature of 30 K. The $\mathrm{H_{2}O}$ o/p ratio has been measured towards about half of them. The similar $\mathrm{H_{2}O}$ and $\mathrm{NH_{3}}$ $T_{\mathrm{spin}}$ suggests that they share a common chemical history and/or they follow a similar thermalization mechanism, while $\mathrm{H_{2}CO}$ either forms or thermalizes differently. Observations of $\mathrm{H_{2}CO}$ o/p ratios in more comets are needed to confirm this is a general behavior.

Towards disks, only the water o/p ratio has been measured so far and only in one source, TW Hya, but the range of estimated ratios is too wide to provide any constraints. Hogerheijde et al. (2011) first reported an $\mathrm{H_{2}O}$ o/p ratio of $0.77\pm 0.07$, but later Salinas et al. (2016) explored different disk models and found values between 0.2 and 3. The disk-averaged $\mathrm{H_{2}CO}$ o/p ratio measured in the HD 163296 disk is consistent with the $\mathrm{H_{2}O}$ o/p ratio measured in TW Hya.

It has long been thought that the o/p ratio upon formation is preserved and $T_{\mathrm{spin}}$ should reflect the formation temperature. This is because the ortho-para inter-conversion by spontaneous radiative transitions and non reactive collisions are extremely slow (e.g. Pachucki & Komasa, 2008; Tudorie et al., 2006), so the o/p ratio upon formation should be preserved.

In general, the o/p ratio upon formation of any molecule will depend on its formation mechanism, the spin selection rules, and the o/p of the reactants. Then, the ortho and para symmetry species can relax to the lowest energy rotational levels through spontaneous emission. Ortho-para conversion can occur via reactive collisions, for example with hydrogen atoms, in which a proton is exchanged. This process, known as equilibration, will compete with the destruction of the molecule. Whether equilibration or destruction dominates will depend on the specific rates, and will determine the resulting o/p ratio of the molecule. A good description of the different interstellar processes that determine the o/p ratio is given by Herbst (2015), for the example of the water cation.

In the specific case of $\mathrm{H_{2}CO}$, the dominant gas-phase formation pathway are reactions of atomic oxygen with CH₃. The destruction of $\mathrm{H_{2}CO}$ will depend on the conditions of the gas. In the presence of strong UV fields (i.e., in the disk surface layers and inner disk), $\mathrm{H_{2}CO}$ will mainly be photo-dissociated, while in UV-shielded regions (i.e., in the disk midplane and outer disk) it will mainly be destroyed by reactions with ions (Guzmán et al., 2011). The original $\mathrm{H_{2}CO}$ o/p ratio will thus depend on the o/p ratio of CH₃. CH₃ forms from the dissociative recombination of CH${}^{+}_{5}$ with electrons, where CH${}^{+}_{5}$ in turn forms from a series of reactions involving H₂, starting from C⁺. The fraction of ortho and para formed in each of these reactions (the branching ratios) will depend on the details of the formation mechanism, which are not fully understood. When a di-hydrogenated molecule like $\mathrm{H_{2}CO}$ forms in the gas-phase through full proton scrambling mechanisms, the resulting o/p ratio should reflect the o/p ratio of its parent molecule. Gas-phase formation in cold gas that is enriched in p-$\mathrm{H_{2}}$ should thus result in o/p ratio $<3$, while formation in the warmer gas should result in the high-temperature limit. However, if full proton scrambling mechanisms are not possible/allowed and only restrictive mechanisms pertain (such as H-abstraction), then even a cold gas-phase formation pathway could result in a statistical ratio of 3 resulting in a complete loss of chemical history information. The exact calculations require an accurate treatment of the selection rules (see for example Gerlich et al. (2006) for the case of H${}_{3}^{+}$, and Le Gal et al. (2017) for the case of H₂Cl⁺). These do not yet exist, and theoretical studies focused on $\mathrm{H_{2}CO}$ are needed to determine which of these processes dominate, and thus whether the thermal history of the gas can be traced with $\mathrm{H_{2}CO}$ o/p ratio measurements.

The resulting o/p ratio from grain surface formation is even less clear, although it is often suggested that the $T_{\mathrm{spin}}$ will reflect the temperature of the ice (e.g., Kahane et al., 1984). Only recently have the nuclear-spin conversion on ices, by thermal and non-thermal desorption processes, been investigated. Laboratory experiments have shown that after formation on water ices (at $<10$ K), $\mathrm{H_{2}}$ has the statistical o/p ratio of 3. But the o/p ratio was found to decrease substantially if $\mathrm{H_{2}}$ is re-trapped in the ice (Watanabe et al., 2010). The resulting $\mathrm{H_{2}}$ o/p ratio thus seems to depend on how long the molecule resides in the ice. More recently, Hama et al. (2016) investigated the o/p conversion of water on ices, and found that the statistical value of 3 is preserved after the water ices are thermally desorbed into the gas-phase (when heating the ices to $\sim 150$ K), no matter if the $\mathrm{H_{2}O}$ ice is formed by vapor-deposition or it is formed in-situ at 10 K. One could argue that the water o/p ratio will equilibrate to the (high) temperature during thermal desorption, while it can remain low in the case of photo-desorption. However, Hama et al. (2016) also investigated the case of photo-desorption, and found that the desorbed water also has an o/p ratio of 3. This unexpected result raised the question of what is the real meaning of the low o/p ratios measured in comets, star and planet forming regions.

Assuming that $\mathrm{H_{2}CO}$ behaves similarly to water, it is possible that the $\mathrm{H_{2}CO}$ o/p ratio is not indicative of the ice formation temperature, but could instead provide information on the desorption mechanism and/or the physical conditions of the gas at the time of desorption. Indeed, the spin temperature inferred from $\mathrm{H_{2}CO}$ in HD 163296 agrees well with the expected gas temperature in the regions where the bulk of $\mathrm{H_{2}CO}$ is expected to form. More laboratory experiments and theoretical predictions are needed to better understand the different formation and desorption processes and how they affect the resulting o/p ratio of $\mathrm{H_{2}CO}$.

Finally, the current SMA observations are sensitive to the colder outer disk, where $\mathrm{H_{2}CO}$ is formed predominantly by grain surface chemistry. Future observations at higher-angular resolution will allow us to measure variations of the $\mathrm{H_{2}CO}$ o/p ratio in the disk. A direct comparison of the o/p ratio with the gas temperature structure of the disk at different radii will help us test different formation scenarios and elucidate what information can be extracted from o/p measurements.