Torsional-rotational spectrum of doubly-deuterated dimethyl ether (CH$_3$OCHD$_2$) -- First ALMA detection in the interstellar medium

In 2013, we have published the first rotational analysis and detection of mono-deuterated dimethyl ether in the solar-type protostar IRAS 16293-2422 with the IRAM 30m telescope. Dimethyl ether is one of the most abundant complex organic molecules (COMs) in star-forming regions and their D-to-H (D/H) ratios is important to understand its chemistry and trace the source history. We present the first analysis of doubly-deuterated dimethyl ether (methoxy-d2-methane, 1,1-dideuteromethylether) in its ground-vibrational state, based on an effective Hamiltonian for an asymmetric rotor molecules with internal rotors. The analysis covers the frequency range 0.15-1.5THz. The laboratory rotational spectrum of this species was measured between 150 and 1500 GHz with the Lille's submillimeter spectrometer. For the astronomical detection, we used the Atacama Large Millimeter/submillimeter Array (ALMA) observations from the Protostellar Interferometric Line Survey, PILS. New sets of spectroscopic parameters have been determined by a least squares fit with the ERHAM code for both symmetric and asymmetric conformers. As for the mono-deuterated species, these parameters have permitted the first identification in space of both conformers of a doubly-deuterated dimethyl ether via detection near the B component of the Class 0 protostar IRAS 16293-2422.


Introduction
By may 2021, around 220 molecules have been detected in the interstellar medium (ISM) or circumstellar shells 1 . Complex organic molecules (COMs) are ubiquitous in the ISM, especially close to the forming low-mass stars where dust temperature is high enough to sublimate ice. With its 9 atoms, dimethyl ether (CH 3 -O-CH 3 , hereafter DME), first detected in the ISM by Snyder et al. (1974) is a large COM of relevance for astrochemistry and one of the main COMs present in these warm and dense inner regions of the envelopes of Class 0 protostars, so called hot corinos (Ceccarelli et al. 2004). With the high sensitivity offered by new observations from the Atacama Large Millimeter/submillimeter Array (ALMA) new opportunities are opened-up for characterising the abundances of these COMs systematically (see, e.g., Jørgensen et al. 2020, for a recent review).
One of the characteristics of the molecules present in the warm gas in hot corinos are their high degrees of deuteration. Although, the cosmic deuterium abundance relative to hydrogen is 1.5−2.0×10 −5 (e.g., Linsky 2003;Prodanović et al. 2010), the abundances of deuterated COMs relative to their non-deuterated counterparts are found to be much higher up to ≈ 10% or more in some cases (e.g., Parise et al. 2006;Jørgensen et al. 2018). 1 https://cdms.astro.uni-koeln.de/classic/molecules These enhancements are thought to be a result of the exothermic reaction H + 3 + HD H 2 D + + H 2 + ∆E that in gas at low temperatures enhances H 2 D + relative to H + 3 . This effect is even more pronounced when one considers the multiple-deuterated variants that for several species (e.g., D 2 CO, Persson et al. 2018; CHD 2 OHCHO, Manigand et al. 2019, D 2 O, Jensen et al. 2021 are even further enhanced compared to the relative abundances of the mono-and non-deuterated species, possibly related to the build-up of the ice-mantles where these species are formed. Dimethyl ether is an interesting target species to extend these studies due to its structure and high abundance. The detection of the mono-deuterated form of CH 3 OCH 3 toward the Class 0 protostellar system IRAS 16293-2422 was first reported by Richard et al. (2013) based on the detection of 20 lines of its symmetric and asymmetric conformers. Through ALMA observations Jørgensen et al. (2018) determined the D/H ratio for these conformers (corrected for the number of equivalent hydrogen atoms) of about 3% toward one component in this source, IRAS 16293B. With the high overall column density of DME toward this source, its D/H ratio and potential enhancement of the doubly deuterated species, this would be a natural place to test new spectroscopic predictions for the multi-deuterated variants.
As we already have detailed in the analysis of the monodeuterated species, DME is a near-prolate asymmetric top (Ray's  Fig. 1: Representation of the doubly-deuterated DME with hydrogen atoms that are pictured in white and deuterium in blue. Figures 1a and 1b show the two possible configurations for the asymmetric conformer while Fig. 1c illustrates the symmetric one. asymmetry parameter κ = −0.922) with only a b-dipole moment component, µ b = 1.302 D (Blukis et al. 1963). The partial deuteration of one CH 3 group leads to two possible configurations. The Figure 1 illustrates the molecule studied in this paper, the doubly-deuterated DME (2D-DME), also called methoxyd 2 -methane in order to avoid confusion on the D atoms locations. The asymmetric conformation is identified when one of the deuterium is located on the C-O-C plane and the second in one of the other two locations. Therefore this conformer has two equivalent configurations (a) and (b) in Fig. 1 with possible tunneling effect between them as detected for some lines of the mono-deuterated species (Richard et al. 2013). When both deuterium atoms are outside of the C-O-C plane, the conformation is called symmetric, as it has a symmetry plane and thus belongs to C s symmetry point group. In this study, the two conformations of 2D-DME were considered as two independent asymmetric top molecules each with a single unsubstituted CH 3 internal rotor. The independence in the treatment of their rotational spectra consisted in ignoring the possible effects of tunneling between two equivalent configurations of the asymmetric conformation or between symmetric and asymmetric conformations. This paper presents the first analysis of the 2D-DME from 0.15 to 1.5 THz (or 150 to 1500 GHz) and reports its detection toward the B component toward IRAS 16293-2422.
Potassium methoxide and p-toluenesulfonyl chloride were purchased from Sigma-Aldrich. Methan-d 2 -ol was purchased from Eurisotop. 2D-DME DME was synthesized by the reaction of potassium methoxide with 4-methylbenzenesulfonic acid, methyl-d 2 ester. This latter was prepared as previously reported for the trideutero derivative but using methan-d 2 -ol as alcohol (Yamamoto et al. 2016). In a three-necked flask under nitrogen and connected to a trap immersed in a cold bath of dry ice were introduced potassium methoxide (5.6 g, 80 mmol) and dry DMSO (30 mL). The mixture was heated to 60°C and 4-methylbenzenesulfonic acid, methyl-d 2 ester (8.2 g, 40 mmol) in dry DMSO (20 mL) was added dropwise. The doubly-deuterated DME was distilled off as it formed and condensed in the cold trap. At the end of the addition, the mixture was heated for 1 hour at 80°C. The trap was then closed with stopcocks. Yield: 1.0 g, 53%.

Lille THz spectrometer
The absorption spectra were measured between 150 and 1500 GHz using the Lille spectrometer Zakharenko et al. (2015). The absorption cell was a stainless steel tube (6 cm in diameter, 220 cm in length). The measurements were performed at typical pressures of 35-60 Pa at room temperature. The different frequency ranges were covered with various active and passive frequency multipliers with the Agilent synthesizer (12.5-18.6 GHz) used as the primary signal source. To increase the sensitivity of the spectrometer, frequency modulation at 20.5 kHz of the reference source and lock-in detection are used. The demodulation of the detected signal may be performed either at 1 f or 2 f , but 2 f demodulation is preferred because of simpler presentation of observed spectrum in this case. Absorption signals were detected by an InSb liquid He-cooled bolometer (QMC Instruments Ltd.). Estimated uncertainties for measured line frequencies are 30 kHz, 50 kHz, and 100 kHz depending on the observed S/N ratio and the frequency range.

Spectral analysis
Similarly to the mono-deuterated species, the spectrum of the 2D-DME shows a very dense structure which is caused by the lines of the ground and the lowest excited vibrational states of symmetric and asymmetric conformations. In addition, due to internal rotation of the unsubstituted methyl top, each rotational level of the two conformations is split into A and E symmetry components of C 3v symmetry point group. Despite relatively high barrier to internal rotation V 3 ≈ 900 cm −1 (see Table 1), the A − E splittings of rotational transitions were resolved for the majority of the ground state lines of the 2D-DME adding thus more complexity to the measured spectra. To treat the internal rotation, we used the ERHAM code (Groner 1997(Groner , 2012. The model used in the code refers to the so-called "combined axes methods". It combines the rho-axis system in which the molecular Hamiltonian is set up, and the principal axes system to which the Hamiltonian is transformed. In the rho-axis system, the molecular internal z axis is set parallel to the ρ vector whose coordinates are calculated using the following expression: where λ g are the direction cosines of the internal rotation axis of the top in the principal axis system, I g are the principal inertia moments and I α is the inertia moment of the methyl top. The magnitude of the rho vector represents the coupling between internal rotation of the methyl top and overall rotation of the molecular frame. For the two conformations of 2D-DME, the rho values are relatively high (see Table 1). They represent an intermediate case of the torsion-rotation coupling in the line between cis methyl formate, HCOOCH 3 , ρ = 0.08 (Ilyushin et al. 2009), and acetaldehyde, CH 3 C(O)H, ρ = 0.3 (Smirnov et al. 2014). For the latter two, strong torsion-rotation coupling, and medium height barriers to internal rotation significantly complicated the analysis and required the inclusion of many higher order terms into Hamiltonian. In the present analysis of the 2D-DME, high barrier to internal rotation simplified the treatment, and in addition to pure rotational Hamiltonian only 5 terms were needed to describe torsional-rotational interaction at close to experimental accuracy.
The symmetric conformation of 2D-DME has a symmetry plane and thus belongs to the C s symmetry point group. Therefore, in the ERHAM code, the symmetry parameters ISCD is set to −1 and α, the angle between the rho vector and ab principal axes plane, is set to 0. On the other hand, the asymmetric conformation has a minimal symmetry and belongs to the C 1 group. Consequently, for the asymmetric conformation, ISCD=+1 and after various attempts to set a value to α, we have decided to fix the parameter to 0 for the same reasons as discussed in the Section 4 of the mono-deuterated analysis (Richard et al. 2013).
The analysis was started with a pattern recognition of a few K a = 0 and K a = 1 a R 0,1 lines in order to produce a fit and a first prediction. Then, the fit was iteratively improved, releasing more parameters, by adding new identified lines.
One difficulty encountered during the analysis was caused by several series of high K a b Q 1,−1 transitions. An example is given in Fig. 3 for K a = 24 series. In this case, the spectral lines in the band head represent a congestion of many rotational transi- tions with different J quantum numbers and different torsional symmetries. Such congestion has a strong influence on spectral line shape. Taking the spectral congestion and the second derivative spectrum into account, it leads to significantly reduced line frequency measurement accuracy, typically 3 to 10 times less than usual. In the weighted least-squares fit implemented in the ERHAM code, the weight of each frequency is calculated as reciprocal of the measurement accuracy squared. The influence of such highly congested lines would be 10 to 100 times less important. Therefore, we decided to remove these congested lines from the final fit and to include only the lines from such series that could be considered as isolated ones.  A total of 1451 lines with J max = 71, K a,max = 27 and 2080 lines with J max = 72, K a,max = 25 were assigned respectively for symmetric and asymmetric conformers. For each conformer, a set of 15 pure rotational and 5 torsional-rotation parameters of the effective Hamiltonian implemented in the ERHAM code was needed to fit the observed molecular lines. The determined parameters are listed in Table 1. In the torsional-rotational part of the Hamiltonian, besides ρ, we determined the following parameters: β the angle between rho vector and principal axis a, ε 10 the tunneling parameter representing the splitting between A and E symmetry sublevels, and two rotational constant tunneling parameters ([A − (B + C)/2] q=1 and [(B − C)/4] q=1 ).

Discussion
This study represents the first comprehensive characterization of the rotational spectrum of the two conformations of 2D-DME (CH 3 OCHD 2 ) in the ground vibrational state, and in the fre-quency range up to 1.5 THz. Experimental measurements, partially given in appendices in Table A.1 and A.2 are available in their entirety in electronic form at the Centre de Données astronomiques de Strasbourg (CDS).
In the ERHAM code, the barrier to internal rotation is not calculated explicitly. Consequently, we used the program BAR-RIER which is derived from the program ASTOR described in Groner et al. (1986) that determines the potential barrier for a single simple internal rotor from torsional transitions or splittings. The ground state energy splitting between the A and E symmetry torsional substates, ∆E(E − A), was combined with the derived value of the internal rotation constant F to obtain an estimate of the potential barrier V 3 for both conformers. These values are direct output of ERHAM. In particular, the internal rotation constants F, is obtained by ERHAM from ρ, β and the rotational constants. The F values determined in this study are 6.223 and 6.338 cm −1 , respectively, for symmetric and asymmetric conformers. The determined V 3 values are given in Ta-ble 1. The estimated error of the V 3 barriers is derived with the assumption that the relative error of V 3 is the same as the relative errors of the ∆E(E − A) energy difference and the tunneling coefficient ε 10 . We found respectively an uncertainty of 3 and 1.5 cm −1 for the symmetric and asymmetric conformers. The determined barrier heights agree well with the results of the study of monodeuterated DME (Richard et al. 2013), and with earlier studies of parent DME (Durig et al. 1976;Lovas et al. 1979).
Under Doppler limited spectral resolution, we didn't observe any additional splittings in the ground state rotational lines of the asymmetric conformation as it could be expected owing to possible tunneling between its two equivalent configurations, and as it was previously observed for monodeuterated DME (Richard et al. 2013). It should be noted that for the latter, the tunneling splittings were barely resolved for a limited number of lines. The lack of additional tunneling splittings for the asymmetric conformer may be explained by the increased mass of the internal rotor which is composed of two deuteriums and one hydrogen for 2D-DME. Provided that the barrier height does not vary significantly under H to D isotopic substitution, the increased internal rotor mass reduces the tunneling probability and consequently the splitting between the two tunneling substates (compared to monodeuterated DME). In this case, one could also expect smaller and unresolvable splittings under the Doppler limit in the rotational spectra. The absence of resolved additional splittings as well as the weighted rms deviations of the fits close to 1 support our initial strategy to treat the two conformations independently as two asymmetric rotors each having a single methyl top.

Predictions
Thanks to the new set of spectroscopic parameters given in Table 1 and using the ERHAM software in its "ERHAM V16g-R3" version (downloaded on the PROSPE website (Kisiel 2001)), we calculated predictions in the JPL catalog format (Pickett 1991) up to 1.5 THz.
Two short examples are provided in Tables B.1 and B.2 from 150.1 GHz to 153.5 GHz. The complete tables are available through the CDS in the catalog format.
Numerical values of overall partition functions were computed at 9 different temperatures using the ERHAM output and are reported in the Table 2. A good approximation of these values can be obtained from the simple rigid asymmetric rotor approximation (Townes & Schawlow 1975) and is expressed in the following formula: In equation 2, the degeneracy factor g n is equal to 1 for the symmetric conformation and 2 for the asymmetric one, since the latter has two equivalent configurations and each energy level is split into two sublevels. We thus obtain 27.53433T 3/2 for the symmetric conformer and 54.48133T 3/2 for the asymmetric one. These values have also been listed in the Table 2 for comparison purposes. However, only the results from the fit, highlighted in bold, were used to derive column densities.

Astronomical search for CH 3 OCHD 2
As an application of the new spectroscopy we searched for the two conformers of CH 3 OCHD 2 in Atacama Large Millime-ter/submillimeter Array (ALMA) observations from the Protostellar Interferometric Line Survey, PILS (Jørgensen et al. 2016). The PILS program is an unbiased molecular line survey of the protostellar system IRAS 16293-2422 with ALMA (project id: 2013.1.00278.S; PI: Jes Jørgensen). PILS covers a continuous frequency range from 329.1 to 362.9 GHz of the main atmospheric window in ALMA's Band 7. The spectral resolution over this full range is ≈ 0.2 km s −1 and the sensitivity ≈ 5 mJy beam −1 km s −1 with an angular resolution of 0.5 .
The advantage of looking for the deuterated species toward IRAS 16293-2422 is that it shows high D/H ratios for many organic molecules compared to, in particular, high-mass star forming regions otherwise targeted in studies of complex organic chemistry (e.g. Jørgensen et al. 2020, and references therein). For example, it was toward this source that a number of multiple deuterated species were identified for the first time (e.g. Loinard et al. 2000;Parise et al. 2002Parise et al. , 2003Manigand et al. 2019) and recently the PILS survey has provided a comprehensive inventory of the deuterated content of the warm complex organic molecules close to the protostar Persson et al. 2018;Manigand et al. 2020).
To look for CH 3 OCHD 2 we here targeted a position offset by 0.5 to the southwest from the "B" component of the system (IRAS16293B) -also the target of a number of the previous PILS papers. The advantage of searching at that position is that the lines are narrow (≈ 1 km s −1 ) and that the opacity of both the main lines of the dust continuum is limited. This enables easy assignments of individual transitions without problems with line overlaps and absorption.
At the targeted scales, the densities should be high enough for the excitation of the molecules to occur under Local Thermodynamic Equilibrium (LTE) conditions. We consequently used synthetic spectra calculated under this assumption to identify the lines of CH 3 OCHD 2 within the PILS range and using them to constrain its column density. The synthetic spectrum method is preferred over, e.g., classical rotation diagrams, as it immediately makes it possible both to search for good matches in the data and to show if lines are predicted at frequencies where no emission is seen. In particular, to confirm the presence of rare species to demonstrate that the latter is not the case is equally important in finding good assigned matches. Also, this method makes it possible to predict the profiles for partially blended lines and also to look for overlaps with previously identified species.
Under the LTE assumption the line profiles are calculated assuming Gaussian line-profiles with constant LSR velocity offsets and line widths. Also, the extent of the source emission (and consequently beam filling factor) is kept constant. For these parameters we adopt the same values as for the single-deuterated variants and most other species in the PILS range -a LSR velocity of 2.6 km s −1 , a line width (FWHM) of 1 km s −1 and source size of 0.5 (assumed as a Gaussian distribution). The choices of the LSR velocity and line width can be tested directly when comparing to the observed spectra. The assumed source size mainly affects the results when optically thick lines are considered: in the case of optically thin emission and when species that thought to be co-existent in the gas it does not affect the derived relative abundances.
With these choices of fixed parameters what remains to be fitted are the column densities of the species and their excitation temperature. For the latter we adopt a temperature of 125 K of the non-and singly-deuterated variants of dimethylether, again under the assumption that the doubly-deuterated variants are present in the same gas. We then perform χ 2 -fitting to the rea- sonably isolated lines to constrain the column densities of the symmetic and asymmetric conformers separately.
Figs. 4 and 5 show the 24 lines predicted in this manner to be the strongest for the asymmetic and symmetric conformers. For the two species good fits can be seen for about 20 transitions -a number of which are well-isolated with peak intensities of up to about 200 mJy beam −1 . The derived column densities are 6.1 × 10 15 cm −2 for the asymmetric and 2.7 × 10 15 cm −2 for the symmetric conformers. A'priori one would expect a factor 2 difference between the two due to the symmetry that occurs when two deuterium atoms replace two of the three hydrogen atoms in a CH 3 -group of dimethyl ether. The factor 2.3 we derive when fitting the two conformers independently is consequently an additional verification of the assignment of the two conformers.
Taken together the derived column densities imply a ratio between the column densities of the doubly-to-mono deuterated forms, N(CH 3 OCHD 2 )/N(CH 3 OCH 2 D) of 15-20% for the two conformers. These imply D/H ratios that are about a factor of 5 higher than those inferred from the singly-to-non deuterated variants once corrected for the number of equivalent hydrogen atoms (e.g., Appendix B of Manigand et al. 2019). The trend of strong enhancements of the multi-deuterated variants also seen for other COMs thus also holds for DME, and provide additional constraints on its formation during the cold phases of star formation.

Conclusion
We have presented here a complete high-resolution study of the doubly-deuterated dimethylether over a decade in the frequency range 150 -1500 GHz (0.15-1.5 THz) in its two conformations, symmetic and asymmetric. We were able to determine an accurate set of spectroscopic parameters that allowed us to reproduce measurements with a standard deviation better than 80 kHz. Thanks to these frequency predictions, both conformers have been detected near the "B" component of the protostellar system IRAS 16293-2422. The doubly-deuterated species is found to be enhanced above the D/H ratio inferred from the mono-to-singly deuterated variants at a level similar to other organic molecules. These enhancements make CH 3 OCH 3 an interesting target for further spectroscopic studies of other multiple deuterated variants, which in turn may shed light on the formation of this complex organic during the cold pre-and protostellar stages.  Fig. 4: Spectra for the 24 transitions of the symmetric conformer of CH 3 OCHD 2 predicted to be the strongest given the derived excitation temperature and column density. In the panels, the red line shows the fit to the CH 3 OCHD 2 transition while the blue lines show the predictions for of all other species identified as part of the PILS survey (see, e.g., Drozdovskaya et al. 2019;