© The Authors 2023

1. Introduction

Polyatomic metal-containing molecules are detected only toward the circumstellar envelopes of evolved stars. Diatomic molecules with Al, Na, or K, and a halogen atom have been detected in the prototypical carbon star envelope IRC +10216 (Cernicharo & Guélin 1987). These species are formed in the hot inner parts of the envelope, close to the star. Metals such as Na, K, Ca, Cr, and Fe are found to survive in the gas phase as neutral and/or ionized atoms in the outer layers of IRC +10216 (Mauron & Huggins 2010), which points to a very rich metal chemistry in the cool outer envelope. Various triatomic metal-bearing species containing Mg, Na, Al, K, Fe, and Ca have been detected in the external layers of IRC +10216 (Kawaguchi et al. 1993; Turner et al. 1994; Ziurys et al. 1995, 2002; Pulliam et al. 2010; Zack et al. 2011; Cernicharo et al. 2019a; Changala et al. 2022). To date, the only metal found to lead to molecules with more than three atoms in IRC +10216 is magnesium. These species are MgCCH (Agúndez et al. 2014; Cernicharo et al. 2019b), HMgNC (Cabezas et al. 2013), MgC₄H, MgC₃N (Cernicharo et al. 2019b), MgC₆H, MgC₅N (Pardo et al. 2021), HMgCCCN, and NaCCCN (Cabezas et al. 2023). The observed abundances for these species are well reproduced adopting the formation mechanism proposed by Dunbar & Petrie (2002). Briefly, ionized metal atoms (M⁺) react with large polyynes (C_2nH₂) and cyanopolyynes (HC_2n + 1N) to form, through radiative association, the cationic complexes MC_2nH₂⁺ and MHC_2n + 1N⁺, respectively, which then recombine dissociatively with electrons yielding as fragments the ²Σ radicals MC_2nH and MC_2n + 1N.

In this Letter we report the detection of four new molecules in the direction of IRC +10216. On the basis of precise quantum chemical calculations, exploration of all possible candidates, and chemical modeling, we assign them to the cations MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺.

2. Observations

New receivers, built within the Nanocosmos project¹ and installed at the Yebes 40 m radiotelescope, were used for the observations of IRC +10216 (α_J2000 = 9^h47^m57.36^s and δ_J2000 = +13 ° 16′44.4″). These observations belong to the Nanocosmos survey of evolved stars (Pardo et al. 2022). Additional details on these observations are given in Appendix A.

3. Results

The latest Q-band Yebes 40 m data between 31 and 50 GHz of IRC +10216 have a sensitivity of 0.15–0.32 mK in the scale for a spectral resolution of 220 kHz, and show a forest of sub-mK unidentified lines. Line identification in this work was done using the catalogs MADEX (Cernicharo 2012), CDMS (Müller et al. 2005), and JPL (Pickett et al. 1998). By January 2023 the MADEX code contained 6464 spectral entries corresponding to different vibrational states and isotopologs of 1765 molecules.

3.1. Detection of B1449 (MgC₄H⁺) and B1447 (MgC₃N⁺)

Among the strongest IRC +10216 unidentified lines in the Q band, seven appear in perfect harmonic relation from J_u = 11 to J_u = 17. They exhibit large intensities, ∼2–3 mK, and are shown in the left panels of Fig. 1. Further exploration of the IRAM 30 m radio telescope 3 mm data reveals nine additional lines, previously unidentified, following the same harmonic relation up to J_u = 33 (see Fig. A.1). The derived line parameters for this series, named B1449, are provided in Table A.1. The 16 lines can be fitted together to the standard Hamiltonian of a linear rotor (ν = 2 BJ_u − 4 D ) with B = 1448.5994±0.0013 MHz and D = 63.45±1.15 Hz. The standard deviation of the fit is 200 kHz (60 kHz for the Q-band lines alone). We checked that values of B divided by 2–6 do not match the observations as many lines would be missing. Moreover, the lines do not show any evidence of spectroscopic broadening (see Fig. 1) as they can be fitted with the standard profile for an expanding envelope with 14.5 km s⁻¹ terminal velocity (Cernicharo et al. 2000). Hence, the carrier of B1449 has a ground electronic state ¹Σ.

Fig. 1. Observations with the Yebes 40 m telescope toward IRC +10216. The left and right panels show the lines of B1449 and B1447 (MgC4H+ and MgC3N+), respectively. The line parameters are given in Tables A.1 and A.2, respectively. The abscissa corresponds to the rest frequency assuming a local standard of rest velocity of −26.5 km s−1. The ordinate is the antenna temperature corrected for atmospheric and telescope losses in millikelvin. The red lines show the fitted line profiles adopting a terminal velocity of expansion of 14.5 km s−1 (Cernicharo et al. 2000).

A second series of lines appear at frequencies very close to those of B1449. The lines are shown in the right panels of Fig. 1 and their line parameters are given in Table A.2. These lines can be fitted to a linear ¹Σ molecule with B = 1446.9380±0.0098 MHz and D = 91±23 Hz (B1447). The standard deviation of the fit is 140 kHz. These results point to two molecules with a ¹Σ electronic ground state, or with a ²Σ electronic ground state but with a very low value of the spin-rotation constant (γ ≤ 0.2 MHz). A value for γ larger than 0.2 MHz would produce a detectable broadening given the high signal-to-noise ratio lines of B1449 lines shown in Fig. 1. We note that the tight harmonic relations in the frequencies of each series of lines, together with the different line intensities of the two series, rule out that they could both arise from the fine structure components from a single open shell molecule (e.g., a ²Σ molecule with a large γ). A detailed inspection of the QUIJOTE line survey (Cernicharo et al. 2021a, 2022, 2023) shows that the lines of B1449 and B1447 are not detected in the cold starless core TMC-1.

Ab initio calculations (see Appendix B and Table B.2) indicate that HCCCNCC and HCCNCCC, both isomers of HC₅N, could fit B1449 and B1447 rather well. However, HC₄NC is detected in TMC-1, but not in IRC +10216 (Cernicharo et al. 2020a). Hence, it is difficult to argue that the carriers of B1449 and B1447, both detected in IRC +10216 but not in TMC-1, are HCCCNCC or HCCNCCC. Moreover, the centrifugal distortion constants for B1449 and B1447 is a factor of two larger than those of HC₅N, HC₄NC, and C₅N⁻, which have a rotational constant similar to that of our series (see Table B.2). This points toward two metal-bearing species. MgC₃N, AlC₃N, NaC₃N, and MgC₄H have rotational and distortion constants very similar to that of B1449 and B1447 (see Appendix B and Table B.2). MgC₃N and MgC₄H, which are ²Σ molecules with γ ∼ 4.3 MHz, have been identified in IRC +10216 (Cernicharo et al. 2019b). Some of the isomers of these two species, MgCCNC, MgNCCC, MgCNCC, and HMgC₄, could fit the value of the rotational and distortion constants of B1449 and/or B1447 (see Table B.2), but they have a value for γ similar to that of MgC₃N, and hence their rotational transitions appear as well-resolved doublets of constant separation in our observations. The same argument applies to the isotopologs of MgC₃N and MgC₄H. The vibrationally excited states of these species will exhibit a ²Π vibronic character and cannot be associated with B1449 and B1447. Some molecules containing sulfur and silicon could also match our series. However, most of them have a ²Π ground electronic state and are not detected in IRC +10216 (e.g., SiCCCN, SiC₄H, and HCCCCS; see Table B.2). Some cations of these species such as SiC₃N⁺ and SiC₄H⁺ have rotational constants close to that of our series. However, they differ by several MHz (0.5%; see Table B.2), and the calculated distortion constants are too low compared with those of B1449 and B1447. Hence, they can be also discarded. C₄S has been detected toward TMC-1 (Cernicharo et al. 2021b), but not toward IRC +10216. SiC₄ has B = 1533.8 MHz (see Table B.2) and a ¹Σ ground state. Its lines are prominent in the Q-band data of IRC +10216, with intensities in the range 20–30 mK. Hence, ²⁹SiC₄ and ³⁰SiC₄ could roughly fit the observed intensities. However, their rotational constants have been derived in the laboratory by Gordon et al. (2000) and do not match those of B1449 and B1447.

Among the possible metal-bearing species, AlCCCN and NaCCCN have a ¹Σ electronic ground state and have been observed in the laboratory (Cabezas et al. 2014, 2019). Although their rotational constants do not fit B1449 and B1447, some of their isomers, such as AlCNCC and AlNCCC or NaNCCC and NaCNCC, could match our series of lines. However, AlCCCN has not yet been detected in our data. NaCCCN has been recently reported by Cabezas et al. (2023), but its lines are very weak and its possible isomers have rotational constants that do not agree with those of our series (see Table B.2). Molecules such as CaC₄H, CaC₃N, KC₄H, and KC₃N produce rotational constants below 1000 MHz (Cabezas et al. 2019). MgC₄ or MgC could fit B1449 or B1447 (see Table B.2). MgC₄ is predicted to have a cyclic structure for its ground electronic state (Redondo et al. 2003a). Nevertheless, it has a linear ³Π electronic state that could fit our series of lines. We estimated a spin-orbit constant of −14.5 cm⁻¹ and a rotational constant of 1427 MHz (see Appendix B), which do not fit either B1449 or B1447. The electronic ground state of MgC is ²Σ and its γ is estimated to be 3–4 MHz producing well-resolved doublets for each rotational transition, which are not seen in our data. The anions MgC₃N⁻ and MgC₄H⁻ could be present as the corresponding neutrals have been detected in IRC +10216 (Cernicharo et al. 2019b). However, their rotational constants are calculated to be too low to fit B1449 or B1447 (see Table B.2).

Prior to this work, HCO⁺ was the only cation known to be present in IRC +10216 (Lucas & Guélin 1990; Agúndez & Cernicharo 2006; Pulliam et al. 2011). No other cation is predicted to be abundant by chemical models. However, ab initio calculations by Bâldea (2020) and in this work show that MgC₄H⁺ has a ¹Σ electronic ground state, a dipole moment of 13.5 D, and rotational and distortion constants that fit perfectly those of B1449 (see Tables 1 and B.1). Moreover, our ab initio calculations for MgC₃N⁺ show that this species also has a ¹Σ electronic ground state, a very large dipole moment of 18.5 D, and rotational and distortion constants that fit those of B1447 very well (see Tables 1 and B.1). The expected error in the calculated rotational constants is lower than 0.1% (see Appendix B). Finally, our ab initio calculations indicate that the rotational constant of MgC₄H⁺ is always slightly larger than that of MgC₃N⁺ independently of the level of theory used in the calculations. Hence, the reverse assignment of B1449 to MgC₃N⁺ and B1447 to MgC₄H⁺ can be ruled out (see Appendix B). Consequently, we conclude that the carrier of B1449 is MgC₄H⁺, while that of B1447 is MgC₃N⁺.

3.2. Detection of B599 (MgC₆H⁺) and B594 (MgC₅N⁺)

In addition to B1449 and B1447, we found 16 lines in harmonic relation with integer quantum numbers from J_u = 26 up to J_u = 41 in the Q-band data of IRC +10216. The lines are shown in Fig. 2. The derived line parameters are given in Table A.3. The derived rotational constants are B = 598.7495±0.0011 MHz and D = 6.13±0.43 Hz (B599). The lines do not show any indication of broadening and we checked that dividing the value of B by 2–4 does not match the observations. Another set of 13 lines were found in harmonic relation with identical quantum numbers and B = 594.3176±0.0026 MHz and D = 4.92±1.16 Hz (B594). The lines are shown in Fig. A.2 and their line parameters are given in Table A.4. The lines of these two series are not observed in TMC-1. Hence, similarly to B1449 and B1447, the carriers could have a ground electronic state ¹Σ, or a ²Σ state with γ ≤ 0.2 MHz. The possible carriers include isomers of HC₇N such as HC₅NC₂, HC₄NC₃, and HC₃NC₄ (see Table B.3). However, HC₆NC (Botschwina et al. 1998) has not yet been detected in IRC +10216, nor in TMC-1. Hence, these species can be discarded as possible carriers of B599 and B594. MgC₆ is also a good potential candidate. However, a cyclic structure is predicted for its ground state (Redondo et al. 2003a). Hence, similar arguments to those used for MgC₄ (see Sect. 3.1 and Appendix B) permit us to discard MgC₆ as the carrier of B599 or B594. After considering all possible carriers, which include isomers of MgC₅N and MgC₆H among other species, we conclude that the best candidate for B599 is MgC₆H⁺, for which we estimate B = 598.2 MHz and D = 6.6 Hz. For B594 we conclude that MgC₅N⁺ is the most likely carrier for which our ab initio calculations provide B = 593.9 MHz and D = 6.6 Hz (see Table B.2). Similarly to the case of B1449 and B1447, the expected error of the calculated rotational constants of MgC₆H⁺ and MgC₅N⁺ is again lower than 0.1% (see Appendix B). The calculated dipole moments of MgC₆H⁺ and MgC₅N⁺ are 18.2 and 23.7 D, respectively (see Table 1).

Fig. 2. Same as Fig. 1, but for the lines of B599 (MgC6H+). The line parameters are given in Table A.3. Green stars under the fitted line profile indicate lines marginally blended with other features but for which the frequency can still be estimated. Red stars indicate lines heavily blended with other features and for which the frequency cannot be estimated.

4. Discussion

The observed spatial distribution of carbon chains observed with ALMA in IRC +10216 consists of a hollow shell with a radius of 15″ (Agúndez et al. 2017). Moreover, Guélin et al. (1993) studied the spatial distribution of MgNC and concluded that this radical also arises from a shell located at 15″ from the star. We thus assume here that the four new Mg-bearing species also arise from the same region. The cusped shape of the lines argue in favor of this. From the observed intensities, and through a rotation diagram analysis, we derived the rotational temperatures and column densities for the four species (see Table 1). We note that for MgC₄H⁺, which is the only species detected in the Q band and in the 3 mm domain, it is possible to fit the Q-band lines (J_u ≤ 17) with a rotational temperature of 11.9±0.7 K, which is similar to that derived for the MgC₃N⁺ lines that cover a similar energy range. However, including the 3 mm lines of MgC₄H⁺ results in a rotational temperature of 17.4±0.3 K, which is similar to the rotational temperatures derived for MgC₆H⁺ and MgC₅N⁺ (see Table 1). Taking into account the high dipole moment of the observed species, it seems very unlikely that collisional excitation alone could reproduce the high rotational temperatures derived from the data. Infrared pumping (Cernicharo et al. 2014; Pardo et al. 2018) may play a role in the population of the high-J lines of these molecules. Infrared pumping is responsible for the detection of HC₅N, HC₇N, and HC₉N in highly vibrationally excited states with energies above 100 K (Pardo et al. 2020, 2023).

The MgC₄H⁺/MgC₆H⁺ and MgC₃N⁺/MgC₅N⁺ abundance ratios are, within a factor of ∼2, similar to those of the corresponding neutral species, MgC₄H/MgC₆H and MgC₃N/MgC₅N (see Table 1). However, the column densities of the cations are 50–80 times lower than those of the corresponding neutral species (see Table 1). This difference is likely related to the formation mechanism of these Mg-bearing species and the higher reactivity of cations compared to neutrals.

To shed light on the formation of these Mg-bearing cations we modeled the chemistry adopting the chemical model presented in Pardo et al. (2021). We updated the initial abundance of Mg to 3 × 10⁻⁶ relative to H to better reproduce the column densities of Mg-bearing molecules observed in IRC +10216. We included the Mg-bearing cations detected here assuming that they are formed by radiative association between Mg⁺ and the corresponding radicals C₄H, C₃N, C₆H, and C₅N, while they are mostly destroyed through dissociative recombination with electrons and reaction with atomic H. It is unlikely that these cations react with H₂ because these processes are calculated to be endothermic (see Appendix C). The most critical parameters are the rate coefficients of radiative association between Mg⁺ and radicals. We therefore performed calculations to compute these rate coefficients for Mg⁺ reacting with C₄H, C₆H, C₃N, and C₅N (see Appendix D).

The abundances calculated with our model are shown in Fig. 3, where it is seen that Mg-bearing cations are indeed expected to be abundant in the outer layers of IRC +10216. The model reflects our finding that the neutral Mg-bearing species are more abundant than their corresponding cations. The column densities predicted for MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺ are 2.1 × 10¹¹ cm⁻², 3.5 × 10⁹ cm⁻², 1.0 × 10¹² cm⁻², and 4.6 × 10¹¹ cm⁻², respectively. The chemical model reproduces the observed column densities within a factor of 4 for the Mg-bearing cations, except for MgC₃N⁺, in which case the calculated column density is about 30 times lower than the observed value (compare with values given in Table 1). This is probably due to uncertainties in the rate coefficients of radiative association calculated for the reactions of Mg⁺ with small molecules such as C₃N which has a low-lying ²Π electronic excited state. It is also interesting to note that the chemical model predicts that cations reach a peak abundance at smaller distances from the central star compared to the corresponding neutrals. High angular resolution observations of these neutral and positively charged Mg–carbon chains will certainly help to characterize the chemistry of magnesium in IRC +10216. Longer Mg-chain cations such as MgC₈H⁺ and MgC₇N⁺ are predicted to have very large dipole moments and reasonable abundances. They could be detected in future improved Yebes 40m data of IRC +10216. Smaller cationic species such as MgCCH⁺ and MgNC⁺ will have very low abundances because the radiative association rates of Mg⁺ with CCH and CN are expected to be very low.

Fig. 3. Calculated abundances for Mg-bearing species, neutral and cationic, resulting from the chemical model described in Sect. 4.

5. Conclusions

We report in this work the detection of four new series of millimeter-wave lines harmonically related in the carbon-rich circumstellar envelope of IRC +10216. From the shape and widths of the lines, which are characteristic of optically thin lines arising in the outer envelope, and from the absence of nearby spectral features that could reveal associated K-components, we conclude that they must arise from four linear molecular species with a ¹Σ electronic ground state.

The sharpness of the U-shaped line edges allows us to derive very accurate values for the spectroscopic constants of these species. The rotation constants of two species, B = 1447 and B = 1449 MHz, are similar to that of cyanodiacetylene, HC₅N, a linear C-chain molecule fairly abundant in IRC +10216. This implies C-chains of similar weight. The distortion constants, D, are twice as large as that of HC₅N, which suggests metal-bearing species. Similarly, the B and D constants of the two remaining species (B = 594 and B = 598 MHz) imply metal-bearing C-chains with a weight similar to cyanotriacetylene, HC₇N.

From a compilation of the spectroscopic constants of two scores of sulfur, silicon, or metal-bearing C-chains, their isomers and their ionic forms, derived from laboratory measurements or from our new or previously published ab initio quantum mechanical calculations, we find a remarkable agreement between B and D of the four new species and those of the cations MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺.

Eight linear Mg-bearing molecules have previously been identified in IRC +10216. They include the neutral species HMgC₃N, MgC₄H, MgC₃N, MgC₆H, and MgC₅N. Chemical modeling of the outer envelope predict that these neutral species form from the dissociative recombination of MgC_nH and MgC_nNH⁺ cations, those cations being themselves produced by the radiative association of polyynes with Mg⁺. The cations MgC_nH⁺ and MgC_nN⁺ are produced by the radiative association of Mg⁺ with the radicals C_nH and C_nN. The abundances derived from the modeling for the neutrals and their cationic species are in good agreement with the observed values, giving confidence that our identifications are correct. This is the first time that cationic metal-bearing molecular species have been detected in space.

Acknowledgments

We thank ERC for funding through grant ERC-2013-Syg-610256-NANOCOSMOS. M.A. thanks MICIU for grant RyC-2014-16277. We also thank Ministerio de Ciencia e Innovación of Spain (MICIU) for funding support through projects PID2019-106110GB-I00, PID2019-107115GB-C21 / AEI/ 10.13039/501100011033, PID2019-106235GB-I00, and PID2021-122549NB-C21.

References

Appendix A: Observations and derived line parameters

The observations in the Q band were carried out with the Yebes 40m radiotelescope. A detailed description of the telescope, receivers, and backends is given by Tercero & López-Pérez (2021). Briefly, the receiver consists of two cold high electron mobility transistor amplifiers covering the 31.0-50.3 GHz range (Q band) with horizontal and vertical polarizations. Receiver temperatures are 16 K at 31 GHz and rise to 25 K at 50 GHz. The backends are 2 × 8 × 2.5 GHz fast Fourier transform spectrometers with a spectral resolution of 38.15 kHz providing the whole coverage of the Q band in both polarizations. For the observations of TMC-1 the original spectral resolution was used while for IRC +10216 the data were smoothed to 229 kHz, which corresponds to a velocity resolution of 1.7 km s⁻¹ at 40 GHz. This spectral resolution is high enough to resolve the U-shaped lines of IRC +10216, which exhibit a full velocity width of 29 km s⁻¹ (Cernicharo et al. 2000). The observations of TMC-1 were performed in frequency-switching mode, while those of IRC +10216 were done in position-switching mode at the same elevation. The main beam efficiency of the telescope varies from 0.6 at 32 GHz to 0.43 at 50 GHz. The Q-band data of TMC-1 and IRC +10216 presented here correspond to 758 h and 696 h, respectively, of on-source telescope time. All data were analyzed using the GILDAS package.²

Further observations at 3 mm presented in this paper were carried out with the IRAM 30m radiotelescope and have been described in detail by Cernicharo et al. (2019a). They correspond to observations acquired during the last 38 years covering the 70–116 GHz range and reaching a confidence level (sigma level) between 0.7 and 3 mK. Examples of these data can be found in Cernicharo et al. (2019a) and Agúndez et al. (2014). The observing mode, in which we wobbled the secondary mirror by ±90” at a rate of 0.5 Hz, ensured flat baselines.

The intensity scale used in this work, antenna temperature (), was calibrated using two absorbers at different temperatures and the atmospheric transmission model ATM (Cernicharo 1985; Pardo et al. 2001). Calibration uncertainties were adopted to be 10 % for both the IRAM 30m data and the Yebes 40m data. Additional uncertainties could arise in the case of IRC +10216 from the line intensity fluctuation with time induced by the stellar infrared flux variability (Cernicharo et al. 2014; Pardo et al. 2018). Pointing corrections were obtained by observing strong nearby quasars and the SiO masers of R Leo. Pointing errors were always within 2-3″.

We used the SHELL method of the CLASS program, within the GILDAS package, which is well adapted for the line profiles observed in circumstellar envelopes. A variable window of 50–100 MHz was used to derive the line parameters in this source, which are given in Tables A.1, A.2, A.3, and A.4 for MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺, respectively. The Q-band lines of MgC₄H⁺ observed with the Yebes 40m telescope are shown in Fig. 1, while those observed at higher frequencies with the IRAM 30m telescope are shown in Fig. A.1. The lines of MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺ are shown in Figures 1 (right panels), 2, and A.2, respectively. In all cases we fixed the terminal expansion velocity for the circumstellar envelope to 14.5 km s⁻¹ (Cernicharo et al. 2000). Most lines of MgC₄H⁺ are free of blending. Two of them (J = 29-28 and J = 31-30) have a marginal overlap with another feature, but a fit to the line profile is still possible. For the other species several lines are partially blended, but it is still possible to derive their frequencies (lines indicated in Figures 2 and A.2 with a green star). Only the J = 34-33 line of MgC₆H⁺ is heavily blended (red star in Figure 2).

Fig. A.1. Observed lines of B1449 (MgC4H+) with the IRAM 30m telescope toward IRC +10216. Line parameters are given in Table A.1. The abscissa corresponds to the rest frequency assuming a local standard of rest velocity of −26.5 km s−1. The ordinate is the antenna temperature corrected for atmospheric and telescope losses in millikelvin. The red lines show the fitted line profiles.

Fig. A.2. Same as Fig. 2, but for the lines of B594 (MgC5N+). The line parameters are given in Table A.4.

Appendix B: Ab initio calculations

We performed ab initio calculations to obtain accurate rotational parameters for all the candidates to be carriers of the series of harmonically related lines found. We also carried out calculations for analogous molecular species for which the rotational parameters are experimentally known. In this manner we were able to scale the calculated values for the target species using experimental-to-theoretical ratios derived from the related species.

The geometry optimization calculations for all the species were done using the coupled cluster method with single, double, and perturbative triple excitations with an explicitly correlated approximation (CCSD(T)-F12; Knizia et al. 2009) and all electrons (valence and core) correlated together with the Dunning’s correlation consistent basis sets with polarized core-valence correlation triple-ζ for explicitly correlated calculations (cc-pCVTZ; Hill et al. 2010). At the optimized geometry, the electric dipole moment was calculated at the same level of theory as that for the geometrical optimization. These calculations were carried out using the Molpro 2020.2 program (Werner et al. 2020). The values for the centrifugal distortion constants were obtained using harmonic vibrational frequency calculations, with the second-order Møller-Plesset perturbation theory method (MP2; Møller & Plesset 1934) and the correlation consistent with polarized valence triple-ζ basis set (cc-pVTZ; Woon & Dunning 1993). These calculations were carried out using the Gaussian16 program (Frisch et al. 2016).

Table B.1 shows the results obtained for MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺, while the dipole moments are shown in Table 1. The rotational constant B was calculated at CCSD(T)-F12/cc-pCVTZ-F12 level of theory, while the constant D was calculated at MP2/cc-pVTZ level of theory. The ab initio values calculated for MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺ are scaled by the experimental-to-calculated ratio of the corresponding parameter of MgC₄H, MgC₃N, MgC₆H, and MgC₅N, respectively.

We followed the same procedure for the HCCCNCC and HCCNCCC species, whose scaled rotational parameters, using HC₅N as reference, are shown in Table B.2. As for MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺, the rotational constants B were calculated at CCSD(T)-F12/cc-pCVTZ-F12 level of theory, while the constants D were calculated at MP2/cc-pVTZ level of theory. Table B.2 also provides the rotational constants of the species discussed in Sect. 3.1.

MgC₄ and MgC₆ are predicted to have ¹A₁ cyclic structures as electronic ground states (Redondo et al. 2003a). Nevertheless, these species also have high energy linear isomers with ³Π electronic ground states. For the ³Π linear isomers of MgC₄ and MgC₆ we carried out the same type of calculations, and chose MgC₄H and MgC₆H as reference species to scale the ab initio parameters (see Tables B.2 & B.3). We estimate a spin-orbit constant A_SO of −14.5 cm⁻¹ for both species. Three series of lines corresponding to the ³Π₂, ³Π₁, and ³Π₀ ladders could be expected, with ³Π₂ having the lowest rotation constant and the lowest energy levels. The rotational constant of linear MgC₄ is ∼22 MHz below those of B1449 and B1447 so that it can be ruled out as a possible carrier of these species. The rotational constant of MgC₆ is predicted to be ∼597 MHz (see Table B.3) and could fit B598 and/or B594. However, we could expect to have other series of lines corresponding to the ³Π₁ and ³Π₀ ladders with decreasing intensity that are not found in our data. The possibility that B594 is the ³Π₂ ladder and the B598 the ³Π₁ ladder can be discarded as the lines of B594 have less intensity than those of B598. Moreover, no other series is found around B = 602 MHz that could be assigned to the ³Π₀ ladder. These arguments, together with the fact that the lowest energy isomers are expected to be cyclic, permit us to discard MgC₄ and MgC₆ as carriers of some of our series of lines.

From a comparison between our calculated and experimental values for species observed in the laboratory, we estimate that the predicted rotational constants for all considered molecules have uncertainties better than 0.1%.

Appendix C: Energetics of the reactions of MgC_nH⁺ and MgC_nN⁺ with H₂

We explored, using quantum chemical calculations, the energetics of all possible channels in the reactions of the cations MgC₄H⁺, MgC₃N⁺, MgC₆H⁺, and MgC₅N⁺ with H₂. The energetics of all the reactions were calculated considering the ZPE of the separated reactants and products at the CCSD/aug-cc-pVTZ level of theory. The results are shown in Table C.1. It can be seen that all the reaction channels are clearly endothermic, which means that these cationic species do not react with H₂ at low temperatures. Moreover, those channels calculated to have a lower endothermicity, such as channel 3 in the reactions MgC₃N⁺ + H₂ and MgC₅N⁺ + H₂, are expected to have quite high barriers because the concerted reactions involve the break of two bonds and the formation of two others. These barriers preclude the reactions from taking place.

Appendix D: Radiative association of Mg⁺ with the radicals C_nH and C_nN

Radiative association in collisions of a metallic ion, M⁺ = Mg⁺, and a polar molecule, A = C_nH or C_nN, can be described by the scheme developed by Herbst (1976), Bass et al. (1981), Herbst (1987), and Klippenstein et al. (1996):

(D.1)

The formation of the initial complex, with rate constant K_F, is dominated by the strong attractive ion-electric dipole interactions, which can be written as

(D.2)

where d is the electric dipole of the molecule, q_ion is the charge of the ion, R is the vector joining the center of mass of A to the ion M⁺ (of norm R), with cos θ = R ⋅ d/Rd. In the present case the C_nH and C_nN molecules are linear with a permanent electric dipole along the molecular axis.

Once the (MA⁺)^* complex is formed, it can either dissociate back, with a rate constant K_B, or emit a photon to estabilize the MA⁺ ion, with a radiative constant K_E. Thus, the total radiative association rate constant, K_RA, is obtained as

(D.3)

K_F(T) is rather fast, especially at low temperature, because of the long-range interactions, and the K_RA rate constant depends strongly on the ratio , which is very sensitive to the molecular properties.

The calculation of the K_F(T),K_B(T), and K_E(T) rates is achieved here using statistical methods similar to those previously used by Bass et al. (1981), Klippenstein et al. (1996), and Dunbar & Petrie (2002). First, the microcanonical rate constants for each individual process are obtained as a function of the collisional energy, E. Second, the thermal rate constants are obtained by averaging over E using a Boltzmann distribution. The method is briefly described here.

According to Klippenstein et al. (1996), the microcanonical rate constant at a given collisional energy, E, and total angular momentum, J, can be expressed as

(D.4)

where X = F, B for formation and back-dissociation, respectively. In this expression, N(E, J) is the cumulative reaction probability (Klippenstein et al. 1996) and ρ_F(E, J) and ρ_B(E, J) are the density of states of M⁺ + A reactants and of the (MA⁺)^* complex, respectively. These densities are rather different and play an important role.

The strong ion-dipole interaction in Eq. (D.2) rapidly reorients the system to a collinear geometry, and hence a simple capture model based on the mono-dimensional radial coordinate can be done, similar to the Langevin model, obtaining a thermal rate constant as

(D.5)

where q_e = 1/4 is the electronic partition function for C₃N, C₄H, and C₅N (for two doublet reactants to form a singlet-state complex, all being Σ states), and q_e = 1/8 for C₆H. In this expression C = q_iond/4πϵ₀, which depends on the electric dipole moment of the molecule A, and μ = m_M⁺m_A/(m_M⁺ + m_A) is the reduced mass.

The back-dissociation microcanonical rate constant takes the simple form

(D.6)

where J_A and ℓ are the angular momenta associated with molecule A and of M⁺ with respect to A (with J = J_A + ℓ), and

(D.7)

The density of states of (MA⁺)^* depends on the binding energy and the frequency of the normal modes, and contributes significantly to slowing down the back-dissociation as it increases.

The radiative rate constant is given by (Klippenstein et al. 1996)

(D.8)

where N = n₁, n₂, ...j, n_M is a collective vibrational quantum number specifying the vibrational states of the n_M normal modes of (MA⁺)^*; P_N(E, J) is the probability distribution of level N; and is the radiative rate of level N corresponding to the sum of the Einstein coefficients over all possible lower states. Due to the high density of states, the sum over N in Eq. (D.8) is performed using a Monte Carlo sampling. In the rigid rotor approximation and summing the Hönl-London factors (Whiting & Nicholls 1974; Whiting et al. 1980) over rotational transitions, the vibrational emission rate no longer depends on J and can be defined as

(D.9)

where the matrix elements in the normal mode approximation can be written as

(D.10)

where only one vibrational quantum transition is considered because the electric dipole component, d_α, is approximated around equilibrium, Q_i = 0, as

(D.11)

The parameters needed to obtain the radiative association rate constant, K_RA, were calculated for MgC_nH⁺ and MgC_nN⁺. The RCCSD(T)-F12 calculations were performed with the MOLPRO package (Werner et al. 2020) using the cvt-f12 basis (Hill et al. 2010). The geometries of all reactants and adducts were optimized, and their normal modes were obtained. The normal modes for MgC₄H⁺(¹Σ) and MgC₃N⁺(¹Σ) are listed in Table D.1. The binding energies obtained, 3.78, 3.54, 3.58, and 3.44 eV (including zero-point energy) for MgC₃N⁺, MgC₄H⁺, MgC₅N⁺, and MgC₆H⁺, respectively, do not vary significantly among the four systems. The dipole moments were calculated through calculations at different electric fields, and the derivatives of Eq. (D.11) were obtained by a fit for each normal coordinate Q_i (listed in Table D.1 for MgC₄H⁺(¹Σ) and MgC₃N⁺(¹Σ)). Finally, the long-range interaction is characterized by the electric dipole moments, 2.89, 2.1, 3.39, and 5.54 debye for C₃N(²Σ), C₄H(²Σ), C₅N(²Σ), and C₆H(²Π), respectively.

The formation rate constants, K_F(T), are very fast due to strong charge-electric dipole interaction (see left panels of Fig. D.1). The differences are proportional to the particular dipole moments of the reactants, which vary at most by a factor of three, which is not appreciable at the scale of the figure.

Fig. D.1. Calculated rates for the radiative association of Mg+ with C3N and C4H. Left panels: Thermal rate constant of formation, KF(T), and radiative association, KRA, for MgC3N+ and MgC4H+, as indicated in each panel. Right panels: Microcanonical rate constants for the back-dissociation, KB(E) = ∑J(2J + 1)*KB(E, J), and the radiative emission averaged in each energy interval, KE(E).

The emission rates, shown in the right panels of Fig. D.1, show a progressive increase in the order MgC₃N⁺, MgC₄H⁺, MgC₅N⁺, and MgC₆H⁺, essentially because the dipole moment of the adduct increases along this series, resulting in an increase of the ∂d_α/∂Q_i derivatives.

The back-dissociation rate, K_B, is dominated by the relatively high density of states, increasing with collision energy, as shown in Fig. D.1. The values of K_B(E) at low energies are very similar (≈ 10⁶ s⁻¹) except for MgC₅N⁺ with K_B(E) = 7 10³ s⁻¹, which presents the lowest normal mode frequencies (of ≈ 54 cm−1), and therefore it has the highest density of states. The case of MgC₆H⁺ is special, showing K_B(E) very similar to the smaller systems of the series studied. The C₆H reactant is a ²Π state, while the MgC₆H⁺ adduct is ¹Σ, and there is a Π state in the vicinity of its equilibrium geometry. The Σ − Π interaction introduces important perturbations in the frequencies of the normal modes.

The K_E(T)/(K_E(T)+K_B(T)) ratio decreases rapidly with collision energy, as the number of open states of the reactants increase, producing an increases in K_B(T). The value at the lowest energy interval is particularly important to determine the radiative association rate constant, where K_E(E) is closer to K_B(E), and the value of the ratio is particularly sensitive to the adduct density of states.

The calculated total radiative association rate constants, K_RA(T), are shown in the left panels of Fig. D.1 and listed in Table D.2 for different temperatures. In general, K_RA(T) increases with the size of the system, being larger for MgC₆H⁺ and MgC₅N⁺ than for MgC₄H⁺ and MgC₃N⁺ according to the trends explained above. For MgC_nN⁺ the rates calculated here are within a factor of 5 of those obtained by Dunbar & Petrie (2002) for the related systems MgHC_nN⁺ for the same n. At 10 K, the rates are 4.5 10⁻¹² and 2.4 10⁻¹¹ cm³ s⁻¹ for MgC₃N⁺ and MgHC₃N⁺, while for n = 5, the rates at 10 K are 5.8 10⁻⁹ and 2.4 10⁻⁸ cm³s⁻¹ for MgC₅N⁺ and MgHC₅N⁺, respectively. These differences are partially due to the higher dipole moments of HC₃N and HC₅N (3.7 and 4.3 debye) compared to those of C₃N and C₅N (2.89 and 3.39 debye). The radiative association rate constant obtained for MgHC_nH⁺ (Dunbar & Petrie 2002) is considerably lower than that obtained here for MgC_nH⁺ (a factor of ≈ 100 for MgHC₄H⁺), simply because HC_nH has no permanent electric dipole, while C_nH does, favoring the formation of the adduct.

All Tables

All Figures

Fig. 1. Observations with the Yebes 40 m telescope toward IRC +10216. The left and right panels show the lines of B1449 and B1447 (MgC4H+ and MgC3N+), respectively. The line parameters are given in Tables A.1 and A.2, respectively. The abscissa corresponds to the rest frequency assuming a local standard of rest velocity of −26.5 km s−1. The ordinate is the antenna temperature corrected for atmospheric and telescope losses in millikelvin. The red lines show the fitted line profiles adopting a terminal velocity of expansion of 14.5 km s−1 (Cernicharo et al. 2000).

In the text

	Fig. 2. Same as Fig. 1, but for the lines of B599 (MgC6H+). The line parameters are given in Table A.3. Green stars under the fitted line profile indicate lines marginally blended with other features but for which the frequency can still be estimated. Red stars indicate lines heavily blended with other features and for which the frequency cannot be estimated.
In the text

	Fig. 3. Calculated abundances for Mg-bearing species, neutral and cationic, resulting from the chemical model described in Sect. 4.
In the text

	Fig. A.1. Observed lines of B1449 (MgC4H+) with the IRAM 30m telescope toward IRC +10216. Line parameters are given in Table A.1. The abscissa corresponds to the rest frequency assuming a local standard of rest velocity of −26.5 km s−1. The ordinate is the antenna temperature corrected for atmospheric and telescope losses in millikelvin. The red lines show the fitted line profiles.
In the text

	Fig. A.2. Same as Fig. 2, but for the lines of B594 (MgC5N+). The line parameters are given in Table A.4.
In the text

	Fig. D.1. Calculated rates for the radiative association of Mg+ with C3N and C4H. Left panels: Thermal rate constant of formation, KF(T), and radiative association, KRA, for MgC3N+ and MgC4H+, as indicated in each panel. Right panels: Microcanonical rate constants for the back-dissociation, KB(E) = ∑J(2J + 1)*KB(E, J), and the radiative emission averaged in each energy interval, KE(E).
In the text