© C. Richard et al. 2021

1 Introduction

By May 2021, around 220 molecules have been detected in the interstellar medium (ISM) or circumstellar shells¹. Complex organic molecules (COMs) are ubiquitous in the ISM, especially close to the forming low-mass stars where the dust temperature is high enough to sublimate ice. With its nine atoms, dimethyl ether (CH₃—O—CH₃, hereafter DME), which was 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, the 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 characterizing 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 is their high degree of deuteration. Although, the cosmic deuterium abundance relative to hydrogen is 1.5 − 2.0 × 10⁻⁵ (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). These enhancements are thought to be a result of the exothermic reaction ${\text{H}}_3^{+}+\text{HD}\leftrightarrows {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2+\Delta E$ that in gas at low temperatures enhances H₂D⁺ relative to H${}_3^{+}$. This effect is even more pronounced when one considers the multiple-deuterated variants, which for several species (e.g., D₂CO, Persson et al. 2018; CHD₂OHCHO, Manigand et al. 2019, D₂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 buildup 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₃OCH₃ 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 mono-deuterated species, DME is a near-prolate asymmetric top (Ray’s 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₃ group leads to two possible configurations. Figure 1 illustrates the molecule studied in this paper, the doubly deuterated DME (2D-DME), also called methoxy-d₂-methane in order to avoid confusion as to the locations of D atoms. The asymmetric conformation is identified when one of the deuteriums is located on the C–O–C plane and the second is located in one of the other two locations. Therefore this conformer has two equivalent configurations (a) and (b) in Fig. 1, with a 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 the 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₃ internal rotor. The independence in the treatment of their rotational spectra consisted in ignoring the possible effects of tunneling between twoequivalent 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.

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.

2 Experimental details

2.1 Preparation of the doubly deuterated DME (methoxy-d₂-methane, 1,1-dideuteromethylether)

Potassium methoxide and p-toluenesulfonyl chloride were purchased from Sigma-Aldrich. Methan-d₂-ol was purchased from Eurisotop.

Doubly deuterated dimethyl ether was synthesized by the reaction of potassium methoxide with 4-methylbenzenesulfonic acid, methyl-d₂ ester. The latter was prepared as previously reported for the trideutero derivative, but using methan-d₂-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, potassium methoxide (5.6 g, 80 mmol) and dry DMSO (30 mL) were introduced. The mixture was heated to 60 °C and 4-methylbenzenesulfonic acid, methyl–d₂ 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 h at 80 °C. The trap was then closed with stopcocks. The yield was 1.0 g, thus 53%.

2.2 Nuclear magnetic resonance (NMR) of the doubly deuterated DME

¹ H NMR (CDCl₃, 400 MHz) δ 3.18 (quint, 1H, ²J_HD = 1.6 Hz, CHD₂); 3.23 (s, 3H, CH₃). ¹³C NMR (CDCl₃, 100 MHz) δ 59.5 (quint, ¹J_CD = 21.6 Hz, ¹J_CH = 139.5 Hz (d), CHD₂); 60.1 (q, ¹J_CH = 140.4 Hz,CH₃).

2.3 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 and 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 were used. The demodulation of the detected signal may be performed either at 1f or 2f, but 2f demodulation is preferred because of the simpler presentation of the 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, 50, and 100 kHz depending on the observed signal-to-noise ratio (S/N) and the frequency range.

3 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 (See Fig 2). 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 the C_3v symmetry point group. Despite a relatively high barrier to internal rotation of V₃ ≈ 900 cm⁻¹ (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, 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: $${\rho }_g=\frac{{\lambda }_g{I}_{\alpha }}{I_g},(g=x,y,z),$$(1)

where λ_g represents the direction cosines of the internal rotation axis of the top in the principal axis system, I_g represents 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₃, ρ = 0.08 (Ilyushin et al. 2009), and acetaldehyde, CH₃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, a high barrier to internal rotation simplified the treatment, and in addition to pure rotational Hamiltonian only five terms were needed to describe a 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 α, 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₁ group. Consequently, for the asymmetric conformation, ISCD = + 1 and after various attempts to set a value to α, we decided to fix the parameter to 0 for the same reasons as discussed in the Sect. 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 ^aR_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 ^bQ_1,−1 transitions. An example is given in Fig. 3 for the K_a = 24 series. In this case, the spectral lines in the band head represent a congestion of many rotational transitions with different J quantum numbers and different torsional symmetries. Such congestion has a strong influence on the spectral line shape. Taking the spectral congestion and the second derivative spectrum into account, it leads to a significantly reduced line frequency measurement accuracy, typically three to ten times less accurate than usual. In the weighted least-squares fit implemented in the ERHAM code, the weight of each frequency is calculated to be reciprocal to the measurement accuracy squared. The influence of such highly congested lines would be ten 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 a series that could be considered as isolated ones.

A total of 1451 lines with J_max = 71 and K_a,max = 27 and 2080 lines with J_max = 72 and K_a,max = 25 were assigned for symmetric and asymmetric conformers, respectively. For each conformer, a set of 15 pure rotational and five 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, ε₁₀ 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}$).

Fig. 2 Spectrum of doubly deuterated DME recorded in Lille. Lower panel: overall spectrum with a frequency range up to 1.5 THz. A zoom to some analyzed blended rotational transitions of the symmetric conformer is shown in the upper right corner.

Fig. 3 Spectrum zoom showing series with K′ = 24. All lines are extremely blended, making the analysis very difficult and uninteresting. Therefore they are not been included in the fit. (a)symmetry number: 0(A), 1(E).

4 Discussion

This study represents the first comprehensive characterization of the rotational spectrum of the two conformations of 2D-DME (CH₃OCHD₂) in the ground vibrational state, and in the frequency range up to 1.5 THz. Experimental measurements, partially given in the appendices in Tables 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 BARRIER 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₃ for both conformers. These values are the 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⁻¹ for symmetric and asymmetric conformers, respectively. The determined V₃ values are given in Table 1. The estimated error of the V₃ barriers was derived with the assumption that the relative error of V₃ is the same as the relative errors of the ΔE(E − A) energy difference and the tunneling coefficient ε₁₀. We found an uncertainty of 3 and 1.5 cm⁻¹ for the symmetric and asymmetric conformers, respectively. 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 did not 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 splittingbetween 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.

5 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 to 153.5 GHz. The complete tables are available in their entirety in electronic form at the CDS.

Numerical values of overall partition functions were computed at nine different temperatures using the ERHAM output and are reportedin 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: $$g_n{\sum }_J(2J+1)e^{-\frac{E_{\text{rot}}(J,K)}{kT}}=g_n\sqrt{\frac{\pi }{ABC}{\left(\frac{k}{h}\right)}^3}{T}^{3/2}.$$(2)

In Eq. (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, which are highlighted in bold, were used to derive column densities.

6 Astronomical search for CH₃OCHD₂

As an application of the new spectroscopy, we searched for the two conformers of CH₃OCHD₂ in 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⁻¹ and the sensitivity ≈5 mJy beam⁻¹ km s⁻¹ 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. 2002, 2003; Manigand 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 (Jørgensen et al. 2018; Persson et al. 2018; Manigand et al. 2020).

To look for CH₃OCHD₂, here, we 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⁻¹) and 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₃OCHD₂ within the PILS range and used them to constrain their column density. The synthetic spectrum method is preferred over classical rotation diagrams, for example, 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, confirming 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 adopted 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⁻¹, a line width (FWHM) of 1 km s⁻¹, 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 are 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 adopted 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 performed χ² -fitting to the reasonably isolated lines to constrain the column densities of the symmetric and asymmetric conformers separately.

Figures 4 and 5 show the 24 lines predicted in this manner to be the strongest for the asymmetric 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⁻¹. The derived column densities are 6.1 × 10¹⁵ cm⁻² for the asymmetric conformers and 2.7 × 10¹⁵ cm⁻² 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₃-group of dimethyl ether. The factor 2.3 we derived 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₃OCHD₂)/N(CH₃OCH₂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 they provide additional constraints on their formation during the cold phases of star formation.

Fig. 4 Spectra for the 24 transitions of the symmetric conformer of CH3OCHD2 predicted to be the strongest given the derived excitation temperature and column density. In the panels, the red line shows the fit to the CH3OCHD2 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; Manigand et al. 2020, and references therein).

Fig. 5 Spectra for the 24 transitions of the asymmetric conformer of CH3OCHD2 predicted tobe the strongest given the derived excitation temperature and column density.

7 Conclusion

Here, we have presented 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, symmetric 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₃OCH₃ an interesting target for further spectroscopic studies of other multiple deuterated variants, which in turn may shed light on the formation ofthis complex organic during the cold pre- and protostellar stages.

Acknowledgements

This work has been supported by the “Programme National Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. This project has received financial support from the CNRS through the MITI interdisciplinary programs. J.C.G. thanks the Centre National d’Etudes Spatiales (CNES) for a grant. The research of J.K.J is supported by the Independent Research Fund Denmark (grant number 0135-00123B).

Appendix A Experimental frequencies

Appendix B Predicted frequencies

References

All Tables

All Figures

	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.
In the text

	Fig. 2 Spectrum of doubly deuterated DME recorded in Lille. Lower panel: overall spectrum with a frequency range up to 1.5 THz. A zoom to some analyzed blended rotational transitions of the symmetric conformer is shown in the upper right corner.
In the text

	Fig. 3 Spectrum zoom showing series with K′ = 24. All lines are extremely blended, making the analysis very difficult and uninteresting. Therefore they are not been included in the fit. (a)symmetry number: 0(A), 1(E).
In the text

Fig. 4 Spectra for the 24 transitions of the symmetric conformer of CH3OCHD2 predicted to be the strongest given the derived excitation temperature and column density. In the panels, the red line shows the fit to the CH3OCHD2 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; Manigand et al. 2020, and references therein).

In the text

	Fig. 5 Spectra for the 24 transitions of the asymmetric conformer of CH3OCHD2 predicted tobe the strongest given the derived excitation temperature and column density.
In the text