© The Authors 2024

1. Introduction

The chemistry of prestellar cores is characterised by high amounts of molecular depletion due to the freeze-out of gas-phase molecules onto interstellar dust grains. At the high densities prevailing in prestellar cores, molecules from the gas phase collide and stick to the grains, building up an ice mantle. Millimetre wave observations of rotational transitions indeed show a molecular abundance decrease for most species towards the core centres (e.g. Willacy et al. 1998; Caselli et al. 1999; Bacmann et al. 2002) where the densities are highest. Additionally, infrared observations have established the presence of several species in the solid phase of the dense interstellar medium (Whittet et al. 1983; Boogert et al. 2015), the details of which are now being revealed by the James Webb Space Telescope (JWST; Yang et al. 2022; McClure et al. 2023; Rocha et al. 2024).

The accretion timescale for CO molecules at a temperature of 10 K, expressed in years, can be estimated as t_CO = 5 × 10⁹/n_H2 where n_H2 is the molecular H₂ density (Bergin & Tafalla 2007). At the densities prevailing in prestellar cores, around 10⁵ − 10⁶ cm⁻³, this timescale is much shorter than typical core lifetimes of around 10⁶ years (Könyves et al. 2015). In the absence of any efficient mechanism returning frozen-out molecules back to the gas-phase, CO and other molecules should be absent from the gas phase, contrary to what is observed in these regions: within 10⁶ years, for a density of 5 × 10⁵ cm⁻³, the molecular abundance should drop by several dozen orders of magnitude, while observations show a much more moderate drop of a few orders of magnitude (Pagani et al. 2012). Therefore, desorption mechanisms allowing the frozen-out molecules to return to the gas phase have to be considered (Leger 1983).

At typical prestellar core temperatures, around 10 K, the available thermal energy is insufficient to enable molecular desorption, but several non-thermal mechanisms have been invoked to account for the observed molecular abundances in the gas phase: UV photodesorption, cosmic-ray sputtering, desorption by cosmic-ray induced UV photons (secondary UV radiation), IR photodesorption, or reactive desorption, in which the exothermicity of grain surface reactions contributes to the release of the reaction’s products in the gas phase. In the last decades, experimental studies have endeavoured to characterise the efficiency of the various desorption processes. Photodesorption rates have been measured under broadband UV irradiation for molecules such as CO, N₂, CO₂, and H₂O by Öberg et al. (2007, 2009b,c) or Muñoz Caro et al. (2010). Further work using monochromatic synchrotron UV radiation provided insights into the mechanisms governing UV photodesorption (e.g. Fayolle et al. 2011; Bertin et al. 2012; Fayolle et al. 2013). These experiments coupled with quadrupole mass spectrometry to detect ejected gas-phase species have revealed that some molecules can be fragmented upon UV irradiation so that the photodesorption yield of the intact molecules may represent only a minor component of the total photodesorbed fragments, as is the case for CH₃OH (methanol), HCOOH (formic acid), CH₃OCHO (methyl formate) (Bertin et al. 2016, 2023) and H₂CO (formaldehyde) (Féraud et al. 2019). As an illustration, the photodesorption yield measured by Bertin et al. (2016) for intact CH₃OH is two orders of magnitude smaller than previously derived from experiments using solely infrared spectroscopy ice measurements (Öberg et al. 2009a), showing that it is therefore essential to monitor the products of the desorption process and not only the loss of the molecule in the solid state. Recent experiments of resonant infrared irradiation of CO and CH₃OH ices suggest that IR photons could also induce desorption with a flux integrated efficiency similar to that of UV photons (Santos et al. 2023, and references therein).

Chemical desorption, or reactive desorption, occurs when the excess energy of a grain surface reaction (typically involving two radicals) can be converted into kinetic energy perpendicular to the surface overcoming the binding energy of the product with the grain (Allen & Robinson 1975). However, the efficiency of this process is difficult to evaluate in situ. Recent experiments have shown that while this mechanism can be efficient in some cases, the nature of the substrate (whether composed of bare graphite or covered with water ice) and its porosity (Dulieu et al. 2013), as well as the number of degrees of freedom of the reaction product plays a major role. In particular, Minissale et al. (2016) did not observe chemical desorption of CH₃OH formed from the reaction CH₃O + H, and estimated that the efficiency of chemical desorption from water ice is at most a few percent. Therefore, in the dense interstellar medium, where the composition of ice mantles is dominated by water ice, chemical desorption is expected to have a low efficiency, for large molecules at least.

Cosmic rays are also believed to contribute to molecular desorption (e.g. Leger et al. 1985; Willacy & Williams 1993; Hasegawa & Herbst 1993; Shen et al. 2004) by sputtering ice species to the gas phase. The cosmic-ray induced desorption process consists of two main contributions: the ‘whole-grain’ heating and the ‘thermal spike’ sputtering (erosion) (Bringa & Johnson 2004). While many studies so far have relied on the formalism proposed by Hasegawa & Herbst (1993) for whole-grain heating, cosmic-ray induced (thermal-spike) sputtering can be investigated experimentally by using swift heavy ions as cosmic-ray analogues. Studies of ice film irradiation by heavy ions carried out by Seperuelo Duarte et al. (2009, 2010) have determined sputtering yields for CO₂ and CO ices by modelling the disappearance of the ice constituents from the solid phase. Of interest to account for the observed gas-phase abundances in the interstellar medium, however, is the amount of molecules that return to the solid phase intact per incident ion. Recent heavy-ion ice irradiation experiments coupled with a quadrupole mass spectrometer (QMS) to monitor the sputtered species in the gas phase (Dartois et al. 2019) have measured this parameter, both for pure ices (e.g. CO₂, CH₃OH, CH₃OOCH₃) and for ice mixtures (CH₃OH/CO₂, CH₃OH/H₂O, CH₃OOCH₃/H₂O). For the species in these studies, a large fraction of the molecules are desorbed intact, and the sputtering yield of a minor species in an ice mixture is close to that of the species making up the bulk of the ice (Dartois et al. 2019).

While the above-mentioned experiments provide us with essential quantitative data regarding the desorption processes, they have been included in few chemical models so far, and it is still unclear whether they can account for observed molecular abundances (Wakelam et al. 2021). Moreover, a large body of experimental data is still missing, even regarding major species in the dense interstellar medium.

Formaldehyde (H₂CO) is a ubiquitous species in molecular clouds and prestellar cores (e.g. Dieter 1973; Sandqvist & Bernes 1980; Mundy et al. 1987), with gas-phase column densities close to that of methanol, around 10¹³ − 10¹⁴ cm⁻² in prestellar cores (Bacmann et al. 2003; Bacmann & Faure 2016) and in protostellar envelopes (Parise et al. 2006). In ice mantles, H₂CO has only been tentatively detected towards protostellar envelopes and its abundance is estimated at around 5%, similar to CH₃OH (Boogert et al. 2015; Rocha et al. 2024). Formaldehyde also represents an important intermediate step towards more complex molecules on grain surfaces (Chuang et al. 2016), or towards complex organic molecule precursors like HCO in the gas phase (Bacmann & Faure 2016). In prestellar cores, it is believed to form chiefly via hydrogenation of CO on grain surfaces (Watanabe et al. 2004; Minissale et al. 2016), as an intermediate step to methanol formation. Gas phase formation mechanisms have also been proposed (Guzmán et al. 2011; Ramal-Olmedo et al. 2021), but distinguishing both formation mechanisms requires detailed knowledge of the physical processes (e.g. desorption), which is lacking.

In this paper we present a study of swift ion irradiation of pure H₂CO ices and determine the sputtering yield of the intact molecule. We calculate the sputtering rate over the cosmic-ray spectrum (from 100 eV/amu to 10 GeV/amu) and use it in a simple model to compare the derived H₂CO gas-phase abundance with both H₂CO observations in a prestellar source and the abundances predicted with other non-thermal desorption mechanisms.

The paper is organised as follows: In Sect. 2 we describe the irradiation experiments of thin H₂CO films carried out at the Grand Accélérateur National d’Ions Lourds (GANIL, Caen, France) and in Sect. 3 the IRAM H₂CO observations towards a prestellar core. In Sect. 4 we present the analysis of both the observational and the experimental data and determine the intact sputtering yield of the molecule. In Sect. 5 we discuss the astrophysical implications of our study and conclude in Sect. 6.

2. Experiments

Pure H₂CO ice irradiation experiments were carried out at the GANIL ion accelerator on the IRRSUD beamline¹ with the Irradiation de GLaces d’Intérêt AStrophysique (IGLIAS) setup (see Augé et al. 2018, for a detailed description) between May 14 and 16, 2021. The vacuum chamber was cryocooled down to a temperature of about 10 K, and held at a pressure of ∼5 × 10⁻¹⁰ mbar. Infrared spectra were collected with a Bruker Vertex 70v Fourier transform infrared spectrometer (FTIR), operating with a HgCdTe-detector and a 1 cm⁻¹ spectral resolution. The angle between the spectrometer beam and the sample surface was 12°. Thin H₂CO films with thicknesses of ∼0.53 μm were condensed on a ZnSe window at normal incidence with a needle placed at a distance of around 15 mm. The pressure during injection was typically 10⁻⁷ mbar. The deposition rate was adjusted to get a good optical quality, and the sample thickness was estimated using the band strengths of Bouilloud et al. (2015) and the interference fringes from the infrared spectra. Gas-phase H₂CO was produced from the thermal decomposition of paraformaldehyde (Sigma Aldrich, purity 95%), heated at 100°C in a homemade oven maintained under secondary vacuum (∼10⁻⁷ mbar). The gaseous species sputtered in the chamber were detected with a MKS Microvision 2 QMS. The branching ratios of the species fragments were evaluated during their injection.

Irradiations were conducted with a ⁸⁶Kr¹⁸⁺ beam of 0.86 MeV/u (74 MeV) corresponding to an electronic stopping power of 2830 eV (10¹⁵ molecules cm⁻²)⁻¹ for an H₂CO ice of 0.81 g cm⁻³ (Bouilloud et al. 2015), as computed with the SRIM-2013 code² (Ziegler et al. 2010). The flux was 5 × 10⁸ ions cm⁻² s⁻¹, and the fluence reached at the end of the experiment was 4.36 × 10¹² ions cm⁻². The corresponding doses are computed and compared to typical interstellar doses in Appendix A. Infrared spectra monitoring the ice composition were taken every two minutes starting from the beginning of the irradiation. The QMS continously scanned masses from 0–80 u, with a scan duration of ∼15 s.

Additionally, a control experiment under the same conditions was carried out in which we irradiated a ∼0.81 μm CO₂ (Air Liquide, purity ≥99.998%) ice film with the same beam, corresponding to an electronic stopping power of 3290 eV (10¹⁵ molecules cm⁻²)⁻¹, also computed with SRIM-2013, for a CO₂ density of 1.0 g cm³ (Satorre et al. 2008). For this experiment, the flux varied between 1.1 × 10⁹ and 2.1 × 10⁹ ions cm⁻² s⁻¹ throughout the irradiation and the final fluence was 10¹³ ions cm⁻².

3. Observations

Observations of one prestellar core, L1689B, were carried out between 2013 and 2019 during several campaigns, at the IRAM 30 m telescope, located at Pico Veleta, Spain. L1689B is situated in the Rho Ophiuchi star-forming region, and is characterised by low central temperatures and high central densities (Steinacker et al. 2016). It shows signs of gravitational contraction and appears dynamically evolved (Lee et al. 2001, 2003; Lee & Myers 2011) despite moderate molecular depletion and deuteration (Bacmann et al. 2003). A wealth of oxygen-bearing molecules, including formaldehyde and more complex organics, have been detected in this source (Bacmann et al. 2003, 2012). The observed coordinates are RA: 16^h34^m48.3^s Dec:−24° 38′04″, corresponding to the peak of the millimetre continuum emission in the source. Several frequency setups were used to target rare isotopes of formaldehyde lines at 2 mm and 1 mm with the Eight MIxers Receivers (EMIR) E1 and E2, respectively. We focused on H₂CO rotational transitions with upper energies typically lower than ∼40 K, as levels of higher energies are unlikely to be sufficiently populated at the low temperatures prevailing in the source (typically 10 K). Because the low-lying transitions of the main isotopologue are optically thick and self-absorbed towards the core centre (Bacmann et al. 2003), we present here observations of the optically thinner isotopologues H₂¹³CO and H₂C¹⁸O, in both ortho and para nuclear-spin symmetries. The receivers were connected to the Fourier transform spectrometer (FTS) at a spectral resolution of ∼50 kHz over a bandwidth of 2 GHz, leading to velocity resolutions of ∼0.10 km s⁻¹ at 2 mm and ∼0.07 km s⁻¹ at 1 mm. The data were taken using the frequency switching procedure with a frequency throw of 7.5 MHz to reduce standing waves. Pointing and focus were typically checked every 1.5 − 2 h on nearby quasars or planets, with additional focus checks carried out shortly after sunset or sunrise. Pointing was found to be accurate to 3 − 4″. Calibrations were performed every 10 to 15 minutes, depending on atmospheric stability. The beam sizes vary from 18″ at 134 GHz to 11″ at 220 GHz. The beam efficiencies were taken following the Ruze formula $B_{\mathrm{eff}}\left(\lambda \right)=B_{\mathrm{eff},0}{e}^{-{\left(\frac{4\pi \sigma }{\lambda }\right)}^2}$, with λ the wavelength, σ = 66 μm the surface rms of the 30 m antenna, and B_eff, 0 = 0.863. The forward efficiencies were 93% and 91% at 2 mm and 1 mm, respectively. The spectra were converted from the antenna temperature scale T_a^* to the main beam temperature scale T_mb using the values of the efficiencies given above.

The observational data were reduced with the Gildas/IRAM software³. The data reduction consisted in folding the spectra to recover the signal from the frequency switching procedure, except for the newer data for which folded spectra are delivered by the system. Then a low-order (typically 3) polynomial baseline was fitted to the spectra over line-less spectral regions and was subtracted before the spectra were co-added. The reduced spectra are presented in Appendix B. Finally, the spectral lines were fitted with Gaussians. The spectroscopic parameters of the observed transitions and results of the Gaussian fits are presented in Table1. The error given in the integrated intensity includes the calibration uncertainty, which was taken to be 15% at 2 mm and 20% at 1 mm.

4. Results

4.1. Experiments

The evolution of H₂CO ice infrared (IR) spectra as a function of fluence is depicted in Fig. 1. At fluence zero, the film is made of pure H₂CO, and neither paraformaldehyde nor CO or CO₂ can be seen in the spectra. As soon as irradiation starts, bands characteristic of CO₂ and CO appear, as well as bands around 2920, 1230, 1100, and 920 cm⁻¹, which we attribute to the formation of an H₂CO polymer (polyoxymethylene, hereafter POM). Conversely, H₂CO bands decrease with fluence. Apart from CO, CO₂, and POM, the infrared analysis shows no evidence of the formation of any other chemical species during the irradiation of H₂CO. The column densities of CO₂, CO, and H₂CO were derived using the band strengths given in Table 2 and the integrated absorbance of the bands at 2343 cm⁻¹, 2139 cm⁻¹, and 1725 cm⁻¹, respectively. The POM column density was calculated using a combination mode at 1105 cm⁻¹ and a band strength of 1.9 × 10⁻¹⁸ cm molecules⁻¹ (Noble et al. 2012). This band strength is not known with precision (Tadokoro et al. 1963; Schutte et al. 1993), and we assumed a factor of 2 uncertainty. An additional uncertainty arises from changes in the refractive index value across strong absorbance bands, resulting in some unphysical negative absorbances after baseline correction (see e.g. the C=O stretch of H₂CO). The derived column density evolution as a function of fluence, with uncertainty zones, is presented in Fig. 2. As can be seen, the POM column density is not accurately determined, but it is consistent with a fast and efficient polymerisation of H₂CO at very low fluence. Thus, at a fluence of 3 × 10¹¹ ions cm⁻¹, the ice is mainly composed of POM.

Fig. 1. IR spectra of irradiated H2CO ice. Left panel: Full range spectra, corrected for baseline, with assigned spectral peaks (see Table 2). Right panel: Focus on the evolution of the four bands used to derive column density of CO2, CO, H2CO, and POM in the ice.

Fig. 2. Column density evolution of CO2, CO, H2CO, and POM ices as a function of fluence. In the case of POM, a factor of 2 uncertainty was assumed for the band strength of the combination mode at 1105 cm−1 (see text), as depicted by the shaded green zone. In the case of H2CO, the band strength of the stretching mode at 1725 cm−1 was also evaluated by Schutte et al. (1993) as 9.6×10−18 cm molecules−1 at 10 K, which corresponds to the highest column density value of the shaded red zone. Otherwise, an uncertainty of 20% was assumed for the band strengths.

In similar experiments on non-polymerisable ices such as CH₃OH and CO₂ (Dartois et al. 2019; Seperuelo Duarte et al. 2009), the column density evolution with fluence can be solved through a simple differential equation involving solely radiolysis and sputtering processes (see e.g. Eq. (3) in Dartois et al. 2019). Here, however, we were not able to solve such a differential equation, due to the additional polymerisation process of H₂CO towards POM, whose exact nature is unknown⁴. In addition, while swift ion irradiation of H₂CO is responsible for the formation of POM, it also destroys POM to CO and CO₂ (see the decrease in POM column density beyond a fluence of ∼3 × 10¹¹ ions cm⁻² in Fig. 2). Therefore, in order to derive the sputtering yield, we resorted to another method based on the determination of the reduction in thin film thickness using the Fabry–Perot interference fringes observed in the infrared spectra (Dartois et al. 2020a, 2023), combined with QMS analysis of the gas composition. We emphasise that in order to validate this method, it was first applied to measure the thickness evolution of an irradiated thin film of CO₂, for which the usual fit of the column density evolution via a differential equation was also possible. As reported in Appendix C, both methods were found to give consistent results, which are also in good agreement with the literature.

In the case of H₂CO, as irradiation proceeds, the ice film becomes a mixture of mainly POM and H₂CO for which the refractive index (n) is not available. Fitting the rigorous expression for the transmission (T) of a thin absorbing film on a thick transparent substrate is therefore not possible. Nevertheless, transmittance spectra follow a sinusoidal pattern, where the period is proportional to the thickness. Thus, we used the following simple expression to fit the transmittance spectra in the ranges of wavenumbers (1/λ) free of molecular absorptions:

$$\begin{array}{c}\hfill T=T_0+Asin\left(\frac{4\pi nd}{\lambda cos({12}^{°})}\right).\end{array}$$(1)

Here T₀ and A are free parameters, n is the refractive index taken as n = 1.33 ± 0.04 (Bouilloud et al. 2015), and d is ice thickness. The factor cos(12° ) corrects for the angle between the infrared beam and the sample. Fits obtained with this expression are shown in Fig. 3. The ice film is composed of POM and H₂CO, with a composition changing during irradiation. Therefore, evaluating the sample mass density ρ is challenging. Furthermore, even the mass density of a pure H₂CO ice film is unknown. Consequently, it is assumed to be constant at the value for liquid formaldehyde, ρ = 0.81 g cm⁻³ (Bouilloud et al. 2015), with an assumed uncertainty of 20%. It was used to derive the surface density of the ice σ = ρd, as a function of fluence, as shown in Fig. 4. The total mass sputtering yield can be directly estimated from the slope of the surface density evolution with fluence as Y_m,tot = (1.8 ± 0.5) × 10⁻¹⁸ g ion⁻¹. It should be noted that this total yield includes all sputtered species from the considered volume, both radiolysed species (C, O, H, H₂, CO, and CO₂) and intact H₂CO fragments.

Fig. 3. Infrared transmittance spectra of a H2CO ice film evolution as a function of fluence. Solid lines: Experimental data. Dashed lines: Model spectra fitted to the data. See text for details.

Fig. 4. Evolution of the surface density as a function of fluence deduced from interference fringes. The fit was performed for fluences greater than 6 × 1011 ions cm−2 since the very beginning of the irradiation corresponds to ice compaction (e.g. Palumbo 2006; Dartois et al. 2013).

In order to determine the chemical composition of the sputtered mass volume, the various fragments ejected into the gas phase were measured with the QMS. The small residual or background signal observed before irradiation was subtracted from the recorded signal during the whole irradiation time, as reported in Fig. 5 (top panel). We note that there is no intact sputtered POM as no signal is observed at m/z = 31, a mass characteristic of POM fragmentation pattern (Butscher et al. 2019; Duvernay et al. 2014). The mass m/z = 30 can be used to follow H₂CO molecules since no fragment at this mass is observed in the CO fragmentation pattern (Fig. 6). Thus, the signals at m/z = 30 and m/z = 28 can be used to retrieve the total number of H₂CO and CO species, respectively, by dividing each m/z signal by its relative contribution in the fragmentation pattern. The remaining signals at m/z = 12 and m/z = 16, after subtracting the contributions of CO and H₂CO, indicate the abundances of C and O atoms, respectively. Once corrected for the total electron-impact ionisation cross-section σ^impact(X) at 70 eV (energy of the QMS electron ionisation source), the relative fraction of each species f(X,g) with fluence F is determined (Fig. 5, bottom panel). QMS measurements at very low masses are not reliable, but the amount of hydrogen (H or H₂) is supposed to be equal to the sum of the amounts of C, O, and CO. We note that CO₂, a minor component of the ice (see Fig. 1), was not detected in the gas phase by the QMS. As irradiation proceeds, the composition of the film changes (Fig. 2) and this is also observed in the gas composition (Fig. 5, bottom panel). The fraction of gaseous H₂CO thus decreases significantly (by approximately 70 %) for initial fluences (up to 2 × 10¹¹ ions cm⁻²), which is consistent with the radiolysis and the rapid polymerisation of solid H₂CO. It then decreases more gradually, before a (noisy) slightly ascending plateau is reached beyond fluences of 3 × 10¹² ions cm⁻². Conversely, the fractions of C and O increase rapidly (by approximately 70%), but also reach a plateau. Overall, the molecular fractions of these three chemical species in the gas vary by less than a factor of two. The larger molecular fractions of CO and H₂ increase up to fluences of about 2 × 10¹² ions cm⁻² before reaching a slightly descending plateau, with variations of less than 15%. We can extract the initial fraction of each species X in the gas phase f_i(X,g) (Table 5), and the subsequent evolution f(X,g) was fitted by third-order polynomials.

Fig. 5. QMS signals as a function of fluence. Top panel: Background corrected signals. Bottom panel: Molecular fractions in the gas (solid lines, left y-axis) derived from the fragmentation patterns and ionisation cross sections, and sputtering yield of each species (dashed lines, right y-axis). See text for details.

Fig. 6. Mass fragmentation pattern monitored during species injection (for ice sample deposition) for H2CO (top panel) and CO (unpublished data) (bottom panel).

By combining the total mass sputtering yield Y_m,tot with the fitted molecular fractions f(X,g), we can derive the sputtering yield of each ejected species Y_X(F),

$$\begin{array}{c}\hfill Y_X(F)=Y_{m,\mathrm{tot}}\times f_m(X,g)\times \frac{N_A}{M(X)},\end{array}$$(2)

where N_A is the Avogadro constant, M(X) is the molecular mass of species X, and f_m(X,g) represents the mass fraction of species X in the gas:

$$\begin{array}{c}\hfill f_m(X,g)=\frac{f(X,g)\times M(X)}{{\displaystyle {\sum }_X}(f(X,g)\times M(X))}.\end{array}$$(3)

The derived sputtering yields for the various ejected species X are plotted as dashed lines in the bottom panel of Fig. 5 (right y-axis). Similarly to the molecular fractions, the intact sputtering yields vary by less than a factor of 2 for C, O, and H₂CO, and by less than 15% for CO and H₂, as expected. The average intact sputtering yields Y_X were computed over the full range of fluences and are listed in Table 5. The average sputtering yield for intact H₂CO is found to be Y_intact = Y_H2CO = (2.5 ± 2.3)×10³ molecules ion⁻¹.

We can now try to estimate the effective sputtering yield, Y_eff. For non-polymerisable species, Y_eff is usually determined using the column density evolution with fluence measured from IR spectroscopy (Dartois et al. 2019; Seperuelo Duarte et al. 2009). For a pure ice X, Y_eff is defined as the semi-infinite thickness effective-sputtering yield for the ice under study. It includes both intact sputtered species X and radiolysed (daughter) sputtered species. Since radiolysis and sputtering occur simultaneously, IR spectroscopy captures the disappearance of H₂CO, regardless of whether it is sputtered intact or radiolysed. Y_eff is simply related to Y_intact as

$$\begin{array}{c}\hfill Y_{\mathrm{eff}}=\frac{Y_{\mathrm{intact}}}{\eta },\end{array}$$(4)

where η is the intact fraction of sputtered H₂CO. In our experiment, the volume to be sputtered just before irradiation is composed of pristine H₂CO ice. The impact of a swift heavy ion induces radiolysis along its track, affecting a fraction of this volume. Therefore, after irradiation has started, the sputtered volume is composed of intact H₂CO, fresh products of radiolysis as well as primary radiolytic products that accumulate in the film (see Fig. 7). Therefore, in order to determine the effective sputtering yield Y_eff, the intact fraction (η) of sputtered H₂CO during irradiation must be estimated using the initial fractions reported in Table 3:

$$\begin{array}{c}\hfill \eta =\frac{f_i({\mathrm{H}}_2\mathrm{CO},g)}{f_i({\mathrm{H}}_2\mathrm{CO},g)+f_i(\mathrm{CO},g)+f_i(\mathrm{C},g)}.\end{array}$$(5)

Fig. 7. Schematic timeline of events during the first impact of a swift heavy ion (SHI) on an ice sample, within the thermal spike model. Chemistry and sputtering occur almost simultaneously within 10−9 s of the ion’s passage (Lang et al. 2020). The length ds corresponds to the sputtering depth probed by an individual ion incident, which is inferior to d0, the initial ice film thickness; rr corresponds to the radiolysis destruction radius, which is at most equal to the sputtering cylinder radius rs (Dartois et al. 2021). Considering the timescales involved, infrared spectroscopy cannot observe the radiolysis effects within the sputtered volume during the passage of an ion.

In practice, we derived an intact fraction of sputtered H₂CO molecules η = 0.23 ± 0.03, corresponding to an effective sputtering yield for the H₂CO film of Y_eff = (10.9 ± 9.9)×10³ molecules ions⁻¹.

In their compendium of selected literature sputtering yields, Dartois et al. (2023) show that the sputtering yield induced by swift heavy ions follows an empirical law that is directly proportional to the square of the electronic stopping power S_e:

$$\begin{array}{c}\hfill Y_{\mathrm{intact}}=Y_0{S}_e^2.\end{array}$$(6)

Here, using a stopping power of 2830 eV (10¹⁵ molecules cm⁻²)⁻¹ for the ⁸⁶Kr¹⁸⁺ beam, we obtain for a pure H₂CO ice Y₀ = (3.1 ± 2.8)×10⁻⁴ 10³⁰ eV⁻² cm⁻⁴ molecule³. However, most studies only have access to the effective sputtering yield Y_eff, which is therefore used for Y_intact in Eq. (6) (e.g. Dartois et al. 2023). In this case, the effective prefactor value is Y₀^eff = (1.4 ± 1.3)×10⁻³ 10³⁰ eV⁻² cm⁻⁴ molecule³.

4.2. H₂CO abundance in L1689B

The H₂CO column density in the L1689B prestellar core was determined using the RADEX radiative transfer code (van der Tak et al. 2007) in the uniform sphere approximation. The collisional rate coefficients used in the computation are the H₂CO−H₂ coefficients from Wiesenfeld & Faure (2013) adapted to the spectroscopy of H₂C¹⁸O and H₂¹³CO and were obtained from the Excitation of Molecules and Atoms for Astrophysics (EMAA) collisional rate database⁵. Because there are no radiative or collisional mechanisms converting one form into the other, the ortho and para nuclear-spin isomers of each isotopologue were treated as two separate species. At thermal equilibrium, the ortho-to-para ratio of H₂CO is ∼1.5 at 10 K, and reaches the statistical value of 3.0 at about 30 K. Only collisions with para-H₂ were considered as the ortho-H₂ abundance represents at most 1 − 10% of the total H₂ (Pagani et al. 2009; Troscompt et al. 2009; Dislaire et al. 2012) and can be neglected in prestellar cores.

Using RADEX, we calculated grids of models for ranges of values in kinetic temperature, H₂ density, and in the molecular column densities (ortho-H₂¹³CO, para-H₂¹³CO, ortho-H₂C¹⁸O, or para-H₂C¹⁸O). The grid points in temperature and H₂ density were chosen to be the same for all species, which is based on the assumption that all of them are coexistent. The densities span a range between 5 × 10⁴ cm⁻² and 5 × 10⁵ cm⁻² and the temperatures between 3 K and 12 K. The column density of ortho-H₂¹³CO N_o13 was varied between 5 × 10¹¹ cm⁻² and 5 × 10¹² cm⁻². The column density of para-H₂¹³CO N_p13 was taken as N_p13 = N_o13/opr where opr is the ortho-to-para ratio in the considered H₂CO isotopologue (either in H₂¹³CO or in H₂C¹⁸O), and the free parameter opr was varied between 1 and 6. For the column density of H₂C¹⁸O, we assumed a ratio cor =  ¹³C/¹⁸O such as the ortho-H₂C¹⁸O column density N_o18 is N_o18 = N_o13/cor and the para-H₂C¹⁸O column density N_o18 is N_p18 = N_o13/(cor × opr), and varied this ratio between 7 and 10 (see e.g. Wilson 1999, for standard values of ¹³C/¹⁸O). The input linewidths for the calculations were taken to be 0.48 km s⁻¹ for H₂¹³CO and 0.38 km s⁻¹ for H₂C¹⁸O, but they have little influence on the results, because our fitting procedure relies on the velocity integrated intensities of the modelled lines compared to that of the observations. The goodness of fit was estimated with a χ² defined as

$$\begin{array}{cc}\hfill {\chi }^2(N_{\mathrm{o}13},T,n_{{\mathrm{H}}_2},\mathit{opr},\mathit{cor})& ={\displaystyle \sum _i}\frac{{(W_i^{\mathrm{mod}}(N_{\mathrm{o}13},T,n_{{\mathrm{H}}_2})-W_i^{\mathrm{obs}})}^2}{{\sigma }_i^2}\hfill \\ \hfill & +{\displaystyle \sum _j}\frac{{(W_j^{\mathrm{mod}}(N_{\mathrm{o}13},T,n_{{\mathrm{H}}_2},\mathit{opr})-W_j^{\mathrm{obs}})}^2}{{\sigma }_j^2}\hfill \\ \hfill & +{\displaystyle \sum _k}\frac{{(W_k^{\mathrm{mod}}(N_{\mathrm{o}13},T,n_{{\mathrm{H}}_2},\mathit{cor})-W_k^{\mathrm{obs}})}^2}{{\sigma }_k^2}\hfill \\ \hfill & +{\displaystyle \sum _l}\frac{{(W_l^{\mathrm{mod}}(N_{\mathrm{o}13},T,n_{{\mathrm{H}}_2},opr,cor)-W_l^{\mathrm{obs}})}^2}{{\sigma }_l^2},\hfill \end{array}$$(7)

where T is the kinetic temperature; n(H₂) is the H₂ density; W^obs is the observed line integrated intensity; σ is the error on the observed integrated intensity; W^mod is the line integrated intensity calculated by RADEX; and the indices i, j, k, l refer to the transitions in ortho-H₂¹³CO, para-H₂¹³CO, ortho-H₂C¹⁸O, and para-H₂C¹⁸O, respectively. Typically, 10 to 15 points were chosen for each parameter. The grids were refined close to parameter values best reproducing the observations. The best fitting parameters are given in Table 4.

Because the ortho-to-para ratio is expected to be at most 3 for H₂CO⁶, we re-ran the model and kept the ortho-to-para ratio set to 3. In this case, the χ² value is slightly worse (5.5 instead of 4.3 for the best model), and the best fitting temperature is slightly lower (closer to 9 K). The other parameters (density, ortho-H₂¹³CO column density, and ¹³C/¹⁸O ratio) remain unchanged. We also ran the fitting routine with the temperature and density values found for CH₃OH by Bacmann & Faure (2016) (i.e. T = 8.5 K and n_H2 = 3.6 × 10⁵ cm⁻³), and again found results for the parameters very similar to the best fit: N_o13 = 1.8 × 10¹² cm⁻².

The line integrated intensities from the best fit model are shown in Fig. 8 against the observed line integrated intensities. The model can account for the observations within the uncertainties.

Fig. 8. Comparison between the observed integrated intensities (blue dots) and the integrated intensities derived from the non-LTE model (red diamonds), for each of the observed lines, labelled by their rest frequencies. In the figure, ‘o’ stands for ortho and ‘p’ stands for para.

With the values in Table 4, we obtain the following values for the molecular column densities: N_o13 = 2.04 × 10¹² cm⁻², N_p13 = 5.4 × 10¹¹ cm⁻², N_o18 = 2.2 × 10¹¹ cm⁻², and N_p18 = 6 × 10¹⁰ cm⁻². Assuming a ratio ¹⁶O/¹⁸O of 557 (Wilson 1999), we obtain a total (ortho + para) column density of 1.6 × 10¹⁴ cm⁻² for the main isotopologue H₂CO.

The H₂ column density at the same position as the H₂CO observations was derived from the model of Steinacker et al. (2016), by integrating the density profile along the line of sight and multiplying by two to account for the contribution of the far side of the core. The uncertainty was determined using the error bars from the fit including the interstellar radiation field (see Steinacker et al. 2016, for details). We find a total H₂ column density N_H2 = 1.3 × 10²³ cm⁻², with a spread 5 × 10²² − 2.1 × 10²³ cm⁻². The corresponding relative abundance of H₂CO with respect to total hydrogen (n_H = n(H)+2n(H₂)) is ∼6 ± 4 × 10⁻¹⁰.

5. Astrophysical applications

5.1. Model

We investigate in this section whether non-thermal desorption from icy grains can explain the gas-phase abundance of H₂CO derived in L1689B. A simple analytical model is employed: H₂CO molecules from the gas phase can only accrete on the icy surface of dust grains, while H₂CO molecules from the grain surface can only desorb via photodesorption (due to both external and internal cosmic-ray induced UV photons) and direct cosmic-ray induced sputtering.

In addition, a steady-state is assumed, with constant kinetic temperature, density, and total H₂CO abundance (gas and ice) throughout the core. This total H₂CO abundance is unknown in L1689B and taken as 3 × 10⁻⁶, which is the typical H₂CO ice abundance in low-mass young stellar objects (LYSOS) (Boogert et al. 2015), corresponding to about 6% with respect to H₂O ice⁷. The kinetic temperature T is fixed at 10 K and the H₂ density at n(H₂) = 2 × 10⁵ cm⁻³, as reported in Table 4 for the L1689B core.

At steady-state, the rates of accretion and desorption are equal and we can write:

$$\begin{array}{c}\hfill R_{\mathrm{acc}}{n}_f^g=R_{\mathrm{des}}{n}_f^s,\end{array}$$(8)

where n_f^g and n_f^s are the formaldehyde densities (in cm⁻³) in the gas and solid ice phase, respectively. The accretion rate (in s⁻¹) can be written as

$$\begin{array}{c}\hfill R_{\mathrm{acc}}=\langle \pi a_g^2\rangle n_g{v}_{\mathit{th}}S,\end{array}$$(9)

where ⟨πa_g²⟩ is the average grain surface (in cm²), n_g the grain number density (cm⁻³), v_th the thermal velocity of H₂CO molecules, and S the sticking coefficient taken as 1.0. The desorption rate (in s⁻¹) can be written as (Flower et al. 2005)

$$\begin{array}{c}\hfill R_{\mathrm{des}}=\frac{{\mathrm{\Gamma }}^{\mathrm{tot}}\langle \pi a_g^2\rangle n_g}{{\displaystyle {\sum }_i}{n}_i^s},\end{array}$$(10)

where ∑_in_i^s is the sum of ice abundances such that n_f^s/∑_in_i^s is the fractional abundance of H₂CO in the ice. In the following, we note that ∑_in_i^s = X_icesn_H. The quantity Γ^tot is the total desorption rate (in molecules cm⁻²s⁻¹) defined as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}^{\mathrm{tot}}={\mathrm{\Gamma }}_{\mathrm{ISRF}}+{\mathrm{\Gamma }}_{\mathrm{CRI}-\mathrm{UV}}+{\mathrm{\Gamma }}_{\mathrm{CRD}},\end{array}$$(11)

where Γ_ISRF, Γ_{CRI − UV}, and Γ_CRD are the desorption rates due to external UV photodesorption, cosmic-ray induced UV phodesorption, and direct cosmic-ray sputtering, respectively. We note that cosmic-ray induced whole-grain heating can be neglected for a thick volatile mantle, as found at high density (n_H ≳ 10⁴ cm⁻³) (Bringa & Johnson 2004; Wakelam et al. 2021).

5.2. Desorption rates

The external UV photodesorption rate can be written as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{ISRF}}=I_{\mathrm{ISRF}}{e}^{-\gamma A_{\mathrm{v}}}{Y}_{\mathit{pd}},\end{array}$$(12)

where I_ISRF is the external interstellar far-UV radiation field (averaged between 6 and 13.6 eV) of ∼10⁸ photons cm⁻² s⁻¹ (Habing 1968); A_v is the visual extinction; γ is a measure of UV extinction relative to visual extinction, which is about 2 for interstellar grains (Roberge et al. 1991); and Y_pd is the photodesorption yield of H₂CO (in molecules photon⁻¹). The last value was recently determined in the laboratory by Féraud et al. (2019), who were able to derive an intact H₂CO photodesorption yield of 4 − 10 × 10⁻⁴ molecule/photon in various astrophysical environments, including interstellar clouds and prestellar cores (see their Table 1). In the following we use an average experimental yield of Y_pd = 7 × 10⁻⁴ molecule/photon for the photodesorption of intact H₂CO.

The cosmic-ray induced UV photodesorption rate can be written as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{CRI}-\mathrm{UV}}=I_{\mathrm{CRI}}{Y}_{\mathit{pd}},\end{array}$$(13)

where I_CRI is now the internal UV field induced by the cosmic-ray excitation of H₂, which was taken as I_CRI ∼ 10³ζ₁₇, where ζ₁₇ is the rate of cosmic-ray ionisation of molecular hydrogen ζ in units of 10⁻¹⁷ s⁻¹ (see the definition of ζ below). This simple relation between I_CRI and ζ₁₇ was deduced from Fig. 4 of Shen et al. (2004), with ζ computed using Eqs. (16)–(17) below.

Finally, the direct cosmic-ray sputtering rate Γ_CRD is computed by summing and integrating the product of the sputtering yield Y_s(ϵ, Z) with the differential cosmic-ray flux j(ϵ, Z):

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{CRD}}=4\pi {\displaystyle \sum _Z}{\displaystyle {\int }_0^{\infty }}2Y_s(ϵ,Z)j(ϵ,Z)dϵ.\end{array}$$(14)

Here ϵ is the kinetic energy per nucleon, Z is the atomic number of the cosmic-ray nuclei, and Y_s(ϵ, Z) is the sputtering yield in H₂CO/ion obtained by combining the experimental yield Y_intact(S_e) for intact sputtering with the calculated electronic stopping power S_e(ϵ, Z). It is multiplied by a factor of 2 to account for the entrance and exit points of the cosmic-rays through the grain. Assuming a quadratic dependence of Y_intact(S_e) with the stopping power, as in Eq. (6), the intact sputtering yield of H₂CO molecules can be described as

$$\begin{array}{c}\hfill Y_s(ϵ,Z)=Y_{\mathrm{intact}}(S_e)=Y_0{S}_e^2=3.1\times {10}^{-4}{S}_e^2,\end{array}$$(15)

where we combined Y_intact = 2.5 × 10³ molecules ion⁻¹ (see Table 3) with the stopping power of 2830 eV (10¹⁵ molecules cm⁻²)⁻¹ for the ⁸⁶Kr¹⁸⁺ beam.

In order to compute Y_s(ϵ, Z) for a range of kinetic energies (see below), the SRIM-2013 code was used to compute the electronic stopping powers S_e for an H₂CO ice density of 0.81 g cm⁻³ and for the cosmic elements with the largest contributions, accounting for the ∼Z⁴ scaling of the sputtering yields (Dartois et al. 2023). In practice, the 15 ions of H, He, C, O, Ne, Mg, Si, S, Ca, Ti, V, Cr, Mn, Fe, and Ni were included using the fractional abundances f(Z) of Kalvāns (2016, their Table 1). It should be noted that iron, despite its low fractional abundance, gives the largest contribution (∼44%) to the rate, as already found by Seperuelo Duarte et al. (2010) for CO ice.

For the differential cosmic-ray flux j(ϵ, Z) (in particles cm⁻² s⁻¹ sr⁻¹ (MeV/amu)⁻¹) we adopted the functional form proposed by Webber & Yushak (1983)

$$\begin{array}{c}\hfill j(ϵ,Z)=\frac{C(Z)E^{0.3}}{{(E+E_0)}^3},\end{array}$$(16)

where C(Z) = 9.42 × 10⁴f(Z) is a normalising constant and E₀ a form parameter that varies between 0 and 940 MeV. This formulation corresponds to a primary spectrum since it neglects the energy-loss of cosmic rays through the dense interstellar gas. It allows, however, the full range of measured ionisation rates to be explored, from standard diffuse clouds (a few 10⁻¹⁶ s⁻¹) to dense clouds (a few 10⁻¹⁷ s⁻¹), by simply varying the parameter E₀. In order to infer the relation between Γ_CRD and ζ, the ionisation rate ζ was computed as in Faure et al. (2019, their Appendix A):

$$\begin{array}{c}\hfill \zeta =4\pi (1+\beta )(1+\mathrm{\Phi }){\displaystyle {\int }_0^{\infty }}{\sigma }_{\mathrm{ion}}^{\mathrm{p}}(ϵ)j(ϵ,1)\mathrm{d}ϵ.\end{array}$$(17)

Here ${\sigma }_{\mathit{ion}}^{\mathrm{p}}(ϵ)$ is the cross-section for ionisation of H₂ by proton impact, Φ is a correction factor accounting for the contribution of secondary electrons to ionisation, and where

$$\begin{array}{c}\hfill \beta ={\displaystyle \sum _{Z\ge 2}}f(Z)Z^2\end{array}$$(18)

is the correction factor for heavy nuclei ionisation. Full details and references can be found in Appendix A of Faure et al. (2019). In practice, the E₀ parameter in Eq. (16) was taken in the range 200−900 MeV so that ζ was explored in the range 8 × 10⁻¹⁸ − 4 × 10⁻¹⁶ s⁻¹. The integrals in Eqs. (14) and (17) were evaluated numerically from ϵ_min = 100 eV/amu to ϵ_max = 10 GeV/amu. The calibration plot between the intact sputtering rate of H₂CO and the ionisation rate is displayed in Fig. 9 (black curve), where Γ_CRD is found to scale linearly with the ionisation rate. We note that although the low-energy content (< 100 MeV) of the cosmic-ray spectrum is poorly constrained, the relation between Γ_CRD and ζ was found to be quite robust.

Fig. 9. Intact sputtering rate of pure H2CO and effective sputtering rates of the pure ices CO, CO2, and H2O as a function of the ζ cosmic-ray ionisation rate. Effective sputtering yield pre-factors are taken from Dartois et al. (2023).

It is now interesting to compare the sputtering rate of H₂CO with those of the main ice components of interstellar grains, that is H₂O, CO, and CO₂, where the latter two have a typical abundance of 20–30% with respect to H₂O. Trace ice species such as NH₃, CH₄, H₂CO, and CH₃OH are expected to be embedded in H₂O-rich, CO-rich, or CO₂-rich layers. We assume that the intact sputtering yield of the trapped molecule is equal to the effective sputtering yield of the ice matrix. Simplistic but intuitive, this hypothesis is actually supported by irradiation experiments performed on methanol diluted in H₂O and CO₂ ice matrices, although the intact sputtering efficiency depends on the dilution factor⁸ (see the detailed discussions in Dartois et al. 2019, 2020b). We note that the effective sputtering yield of either water or CO₂ was similarly employed to model the cosmic-ray induced sputtering of CH₃OH in cold molecular cores (Wakelam et al. 2021; Paulive et al. 2022). The above procedure was thus applied to pure CO, CO₂, and H₂O ices using the effective sputtering yield pre-factors reported in Table 1 of Dartois et al. (2023) with ice densities (needed to compute S_e) of 0.88 g cm⁻³, 1.0 g cm⁻³, and 0.94 g cm⁻³, respectively (Luna et al. 2022; Satorre et al. 2008; Bouilloud et al. 2015). The corresponding sputtering rates as functions of ζ are shown in Fig. 9. As expected, our results for CO, CO₂, and H₂O are in very good agreement with those of Dartois et al. (2023, see their Fig. 5). We can also observe that the intact sputtering rate of pure H₂CO is a factor of 10–100 lower than the effective sputtering rates of H₂O, CO, and CO₂. For convenient use in astrochemical models, simple linear fits of the sputtering rates shown in Fig. 9 were performed over the range ζ = [8 × 10⁻¹⁸ − 10⁻¹⁶] s⁻¹. These fits are reported in Table 5.

5.3. H₂CO gas-phase abundance

Because the total abundance of formaldehyde n_g^tot is conserved, we have

$$\begin{array}{c}\hfill n_f^g+n_f^s=n_f^{\mathrm{tot}}.\end{array}$$(19)

By combining Eqs. (8–10) and (19), we obtain

$$\begin{array}{c}\hfill n_f^g=\frac{n_f^{\mathrm{tot}}}{1+R_{\mathrm{acc}}/R_{\mathrm{des}}}=\frac{n_f^{\mathrm{tot}}}{1+v_{\mathit{th}}{X}_{\mathrm{ices}}{n}_{\mathrm{H}}/{\mathrm{\Gamma }}^{\mathrm{tot}}}.\end{array}$$(20)

Since in practice v_thX_icesn_H/Γ^tot ≫ 1 and substituting $v_{\mathit{th}}=\sqrt{8k_{B}T/(\pi \mu )}$, where μ is the mass of H₂CO, the H₂CO gas-phase fractional abundance, X_f^g = n_f^g/n_H, takes the simple form

$$\begin{array}{c}\hfill X_f^g=\frac{X_f^{\mathrm{tot}}{\mathrm{\Gamma }}^{\mathrm{tot}}}{\sqrt{8k_{B}T/(\pi \mu )}{X}_{\mathrm{ices}}{n}_{\mathrm{H}}},\end{array}$$(21)

where X_f^tot = n_f^tot/n_H is the H₂CO total fractional abundance, taken as 3 × 10⁻⁶ (see above). The total fractional abundance of ices, X_ices, is unknown in L1689B, but again we can assume typical ice abundances in LYSOS, as compiled by Boogert et al. (2015, their Table 2). By including ‘securely’, ‘likely’, and ‘possibly’ identified species, we obtain X_ices ∼ 8 × 10⁻⁵. Finally, using the physical conditions derived from H₂CO emission in L1689B (i.e. T = 10 K and n_H = 4 × 10⁵ cm⁻³), we obtain numerically:

$$\begin{array}{c}\hfill X_f^g=n_f^g/n_{\mathrm{H}}=1.1\times {10}^{-11}{\mathrm{\Gamma }}^{\mathrm{tot}}.\end{array}$$(22)

In order to reproduce the measured abundance in L1689B (∼6 × 10⁻¹⁰), a total desorption rate Γ^tot ∼ 50 molecules cm⁻²s⁻¹ is required. We compare in Fig. 10 the observational abundance to our predictions from the above analytical model, assuming the cosmic-ray ionisation rate is ζ = 3 × 10⁻¹⁷ s⁻¹ (i.e. ζ₁₇ = 3), which is the typical ionisation rate in molecular clouds, where A_v > 10 (see Fig. 19 in Indriolo & McCall 2012). Five different models are presented: desorption is either due to external and internal photons only (in orange) or to both photons and direct cosmic-ray sputtering, assuming H₂CO ices are pure (in black) or embedded in H₂O-rich (in blue), CO₂-rich (in red), or CO-rich (in green) matrices. We can observe that photodesorption alone cannot explain the gas phase abundance of H₂CO in L1689B and that direct cosmic-ray sputtering is much more efficient, except in the case of a pure H₂CO ice. In particular, cosmic-ray sputtering of H₂CO from CO-rich or CO₂-rich ices can reproduce the gas phase abundance measured towards L1689B to within a factor of 3. Interestingly, this result is consistent with observational evidence that catastrophic CO freeze-out occurs at A_v ≥ 9 and n_H > 10⁵ cm⁻³, where CO hydrogenation may form a CH₃OH-rich layer including H₂CO and new CO₂ molecules (Boogert et al. 2015). It is also possible that part of the interstellar H₂CO forms in the gas phase, for example via the barrierless reaction CH₃ + O (see e.g. Terwisscha van Scheltinga et al. 2021). Thus, while the present results emphasise the potential role of direct cosmic-ray sputtering of H₂CO molecules from CO-rich or CO₂-rich ices, the importance of gas-phase formation and destruction of H₂CO will need to be explored in future works. Finally, we note that the H₂CO gas-phase abundance in L1689B corresponds to about 0.02% of the total (gas and ice) H₂CO abundance.

Fig. 10. H2CO gas phase abundance with respect to H nuclei as a function of visual extinction. The observational value derived in this study for the prestellar core L1689B (blue hatched zone) is compared to five different models for H2CO sputtering: pure H2CO ice and H2CO ideally diluted in H2O, CO2, and CO matrices. See text for details.

In summary, our simple model shows that the yields for experimentally determined direct cosmic-ray sputtering are large enough to release significant amounts of H₂CO ice into the gas-phase, provided that H₂CO is embedded in CO-rich or CO₂-rich ice mantles and that H₂CO molecules co-desorb with the matrix (i.e. with the same yield). Obviously, further experiments with mixed ices (e.g. H₂CO:CO), as well as more comprehensive (time-dependent and density-dependent) astrochemical models are needed to confirm these results.

6. Conclusions

We have studied the sputtering of pure formaldehyde ices under cosmic rays using swift-ion irradiation experiments of H₂CO ice films at 10 K carried out at GANIL with the IGLIAS setup. Formaldehyde is found to quickly polymerise upon swift ion exposure. Other detected radiolysis products in the ice are CO and CO₂. We determined the intact H₂CO sputtering yield in the gas phase to be Y_intact = 2.5 × 10³ molecule ion⁻¹ for a stopping power of 2830 eV/10¹⁵ molecule cm⁻². Assuming a quadratic dependence of the sputtering yield with stopping power, this corresponds to an intact sputtering prefactor Y₀ = 3.1 × 10⁻⁴ 10³⁰ eV⁻²cm⁻⁴molecule³. The effective sputtering prefactor is Y₀^eff = 1.4 × 10⁻³ 10³⁰ eV⁻²cm⁻⁴molecule³. With these values we determined that the direct cosmic-ray sputtering rate for intact H₂CO in pure ices can be expressed simply by Γ_CRD(H₂CO) = 0.15 ζ₁₇ molecules cm⁻² s⁻¹, where ζ₁₇ is the cosmic-ray ionisation rate in units of 10⁻¹⁷ s⁻¹.

Using a simple steady-state chemical model including accretion of molecules from the gas phase onto the grains and their desorption by either UV photons (both external and secondary UV photons produced by cosmic-ray excitation of H₂) or direct cosmic-ray impact, we find that photodesorption is not efficient enough to account for the gas-phase abundance of H₂CO of 6 × 10⁻¹⁰ (with respect to H nuclei) measured in the prestellar core L1689B, nor is cosmic-ray induced sputtering of a pure H₂CO ice. On the other hand, when H₂CO is diluted in CO or in CO₂ ices, our model can reproduce the H₂CO gas-phase abundance, assuming that minor ice constituents like formaldehyde will be sputtered to the gas phase at the same effective rate as the ice mantle bulk. However, the cosmic-ray sputtering rate is insufficient if H₂CO is embedded in a water ice mantle. Future experimental studies should be dedicated to investigating the polymerisation and sputtering of H₂CO diluted in other ice matrices (e.g. H₂O, CO, CO₂). The competition between the cosmic-ray induced sputtering and the gas-phase chemistry of H₂CO will also need to be explored.

Data availability

The experimental data can be downloaded from https://zenodo.org/records/13383307

Acknowledgments

This work is based on observations carried out under project numbers 090-17 and 135-18 with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research has made use of spectroscopic and collisional data from the EMAA database (https://emaa.osug.fr and https://dx.doi.org/10.17178/EMAA). EMAA is supported by the Observatoire des Sciences de l’Univers de Grenoble (OSUG). The experiments were performed at the Grand Accélérateur National d’Ions Lourds (GANIL) by means of the CIRIL Interdisciplinary Platform, part of CIMAP laboratory, Caen, France. We thank the staff of CIMAP-CIRIL and GANIL for their invaluable support. We acknowledge funding from the National Research Agency, ANR IGLIAS (grant ANR-13-BS05-0004). Thibault Launois is acknowledged for his contribution to the experiments and data reduction.

References

Appendix A: Comparison of the experimental and interstellar doses

In our experiments, the dose D (D = S_e ⋅ F where F is the fluence) lies between 0.0025 and 12 eV/molecule at the investigated experimental fluences (8.8 × 10⁸ − 4.4 × 10¹² ion cm⁻²). These values can be compared to the "interstellar" dose received by a grain over 1 Myr (the typical age of a molecular cloud or prestellar core), which we estimate as ∼0.3 eV/molecule using the cosmic-ray flux given in Eq. 16, with E₀ = 500 MeV and S_e computed with the SRIM-2013 code. As a result, the most astrophysically relevant conditions in our experiment should be the early times of the measurements, when the fluence is lower than a few ∼10¹¹ ions cm⁻². We note, however, that the experimental and calculated doses correspond to different kinds of ions, Kr and essentially H, He, C, O, and Fe, respectively. It is unclear whether the nature of ions has an effect on radiolysis. Thus, threshold effects have been reported for polymers that contain aromatic species, but not in the case of aliphatic polymers that are chemically more similar to ices (Balanzat et al. 1995). To our knowledge, no similar studies have been published for ices, and firm conclusions about the relevance of experimental doses would require comparative radiolysis investigations with e.g. H, He and Fe ions.

Appendix B: H₂CO spectra

The spectra observed with the IRAM 30 m radiotelescope and used to determine the H₂CO column density are presented in Fig. B.1 for H₂¹³CO and in Fig. B.2 for H₂C¹⁸O. In the figures, ‘o’ and ‘p’ stand for ortho and para spin symmetry, respectively. The spectra were taken at coordinates RA = 16^h34^m48.3^s, Dec = −24° 38′04″.

Fig. B.1. Observed H213CO spectra. The vertical dotted line shows the source velocity (v = 3.5 km s−1). The spin symmetry is indicated in the upper left corner of each plot, with ‘o’ standing for ortho and ‘p’ standing for para.

Fig. B.2. Observed H2C18O spectra. The vertical dotted line shows the source velocity (v = 3.5 km s−1). The spin symmetry is indicated in the upper left corner of each plot, with ‘o’ standing for ortho and ‘p’ standing for para.

Appendix C: Sputtering yield determination: The case of a thin CO₂ ice film

High-energy heavy ions irradiation of CO₂ ices has been thoroughly investigated (Seperuelo Duarte et al. 2009; Mejía et al. 2015; Dartois et al. 2021). In this appendix, we compare sputtering yields estimated from the interference fringes analysis method (Method 1, as employed in this work for H₂CO ice) to yields derived from the integration of IR bands (Method 2).

C.1. Method 1: Interference fringes method

The procedure follows the steps explained in Section 4.1 for H₂CO ice. The interference fringes observed in the IR spectra during irradiation are depicted in Fig. C.1. The adjustments were achieved using Eq. 1. While a more rigorous equation including the optical indices could be applied in the case of CO₂ (see Dartois et al. 2023, for details), we opted for the simplified equation (Eq. 1) in order to be consistent with the methodology employed in this study. The total mass sputtering yield (Y_m, tot = (5.7 ± 0.1)×10⁻¹⁸ g ion⁻¹) was extracted from the linear fit of surface density as a function of fluence (Fig. C.2), using a density of CO₂ of 1.0 g cm⁻³ (Bouilloud et al. 2015). We note that in the case of CO₂, there is no visible compaction phase, despite experimental conditions close to those of H₂CO (similar film thicknesses, pressure and temperature, see Section 2).

Fig. C.1. Infrared transmittance spectra of a CO2 ice film evolution upon 74 MeV 86Kr18+ swift heavy ion irradiation. Solid lines: Experimental data. Dashed lines: Model spectra fitted to the data, without considering the absorption bands. The QMS signal ends around a fluence of 4.1012 ion cm−2, before the end of the irradiation as the instrument stopped prematurely.

Fig. C.2. Evolution of the bulk mass per area as a function of fluence deduced from interference fringes for a CO2 ice film.

Fig. C.3. QMS signals as a function of fluence. Top panel: Background corrected signals. Bottom panel: Molecular fractions in the gas (solid lines, left y-axis) derived from the fragmentation patterns and ionisation cross sections, and sputtering yield of each species (dashed lines, right y-axis). See text for details.

The CO₂ fragmentation pattern was derived from the QMS signal during CO₂ injection. It is composed of only four fragments at m/z = 44, m/z = 28, m/z = 16, and m/z = 12 (with intensity ratios of 0.748, 0.112, 0.07, and 0.07, respectively). Figure C.3 shows the molecular fraction deduced from the QMS signal after (i) correction of the background, (ii) extracting the original CO₂ signal intensity using the fragmentation pattern (I_CO2 = I(m/z = 44)/0.748), (iii) subtracting the CO₂ fragments from other mass signals (e.g. I_C = I(m/z = 12)−0.07 × I_CO2), and (iv) dividing by the total electron-impact ionisation cross sections (σ^impact(CO₂) = 3.521Ų, from the NIST database ).

In contrast to H₂CO, no significant variation of the CO₂ fraction is observed during the very early fluences; however, a similar smooth and continuous evolution is observed thereafter. Over the entire experiment, the molecular fraction of CO increases by ∼40%, while that of CO₂ decreases by ∼60%. The sputtering yields were calculated using Eq. 3 after fitting the molecular fractions with a third-order polynomial, as in the case of H₂CO. The fits are shown as dashed lines in the bottom panel of Fig. C.3. The average sputtering yields are provided in Table 5.

As in the case of H₂CO, the intact fraction of sputtered CO₂ was computed using the initial fractions f_i(X,g) of species containing carbon atom (Table C.1) :

$$\begin{array}{c}\hfill \eta =\frac{f_i({\mathrm{CO}}_2,g)}{f_i({\mathrm{CO}}_2,g)+f_i(\mathrm{CO})+f_i(\mathrm{C},g)}.\end{array}$$(C.1)

We obtain η = 0.25±0.02 from which the effective sputtering yield (see Eq. 4) for a CO₂ ice film was estimated as Y_eff=(6.8±4.1)×10⁴ molecules ions⁻¹. The corresponding effective prefactor derived from Eq. 6 is Y₀^eff=(3.8±2.3)×10⁻³ 10³⁰ eV⁻² cm⁻⁴ molecule³, using a stopping power of 3290 eV/10¹⁵ molecules cm⁻².

C.2. Method 2: Integrated bands method

Fig. C.4 illustrates the evolution of the CO₂ infrared vibrational band at 2343 cm⁻¹ with fluence. The integrated intensity of this band (A) was employed to derive the CO₂ column density (N) with the first-order relation:

$$\begin{array}{c}\hfill N=\frac{ln(10)}{\mathcal{A}}{\displaystyle \int A}(\nu )\mathrm{d}\nu ,\end{array}$$(C.2)

where the band strength (𝒜) is 7.6 × 10⁻¹⁷ cm molecules⁻¹ (Bouilloud et al. 2015). The derivative of the column density with respect to fluence was fitted using an equation involving both radiolysis and sputtering processes:

$$\begin{array}{c}\hfill -\frac{\mathit{dN}}{\mathit{dF}}={\sigma }_{\mathit{des}}N+Y_{\mathit{eff}}(1-e^{-\frac{N}{N_d}})f,\end{array}$$(C.3)

where σ_des is the destruction cross-section of CO₂, f is the relative concentration of carbon dioxide molecules with respect to the total number of molecules in the ice film, and N_d is the semi-infinite sputtering depth. Similarly to the experiment of Mejía et al. (2015) on CO₂ ice irradiation, the fraction of CO₂ in the ice remains significantly higher than that of CO at the observed fluences, and the value of f is set to 1. Adjusting the data with Eq. C.3, as shown in Fig. C.5, provides the following solution: Y_eff = (7.0 ± 0.9)×10⁴ molecules ions⁻¹, σ_des = (1.4 ± 0.1)×10⁻¹³ cm² ion⁻¹, and N_d = (1.8 ± 0.3)×10¹⁷ molecules cm⁻². The derived effective sputtering yield corresponds to an effective prefator Y₀^eff = (6.5 ± 0.8) × 10⁻³ 10³⁰ eV⁻² cm⁻⁴ molecule³, using again a stopping power of 3290 eV/10¹⁵ molecules cm⁻².

C.3. Comparison of Method 1 and 2

Results from this work and previous studies are presented in Table C.2. The comparison with literature data is not straightforward due to various experimental conditions. For example, experiments have shown that the sputtering efficiency depends on the ice thickness: the effective sputtering yields Y_eff for CO and CO₂ ices were indeed found to be constant for thick films (more than ∼ 100 monolayers) but to decrease significantly for thin ice films (Dartois et al. 2021). We recall that Y_eff scales with the square of the heavy ion stopping power, while Y_e⁰ should be roughly constant for a given target. We can observe in Table C.2 that our Y_e⁰ value derived from the integrated band method (Method 2) agrees with those of Seperuelo Duarte et al. (2009) and Dartois et al. (2021), who used the same method, to better than a factor of 2, and almost within error bars (whatever the thickness of the film). In addition, and more importantly, the interference fringes method (Method 1) gives an effective sputtering yield prefactor (Y_e⁰ value) in very good agreement (i.e. within error bars) with the value we derived from Method 2.

Fig. C.4. Evolution of the background-corrected IR spectrum centered on the ν3 stretching mode of CO2 upon 74 MeV 86Kr18+ swift ion irradiation

Fig. C.5. Column density derivative with respect to fluence. The solid and dashed lines are fits with or without including the sputtering depth Nd in Eq. C.3, respectively.

We conclude that the interference fringes method, even using a simplified optical model, is a powerful alternative to the integrated bands method for the determination of sputtering yields. Furthermore, combined with QMS measurements, this method allows us to determine both effective and intact sputtering yields, this latter being challenging to estimate with other methods.

All Tables

All Figures

	Fig. 1. IR spectra of irradiated H2CO ice. Left panel: Full range spectra, corrected for baseline, with assigned spectral peaks (see Table 2). Right panel: Focus on the evolution of the four bands used to derive column density of CO2, CO, H2CO, and POM in the ice.
In the text

Fig. 2. Column density evolution of CO2, CO, H2CO, and POM ices as a function of fluence. In the case of POM, a factor of 2 uncertainty was assumed for the band strength of the combination mode at 1105 cm−1 (see text), as depicted by the shaded green zone. In the case of H2CO, the band strength of the stretching mode at 1725 cm−1 was also evaluated by Schutte et al. (1993) as 9.6×10−18 cm molecules−1 at 10 K, which corresponds to the highest column density value of the shaded red zone. Otherwise, an uncertainty of 20% was assumed for the band strengths.

In the text

	Fig. 3. Infrared transmittance spectra of a H2CO ice film evolution as a function of fluence. Solid lines: Experimental data. Dashed lines: Model spectra fitted to the data. See text for details.
In the text

	Fig. 4. Evolution of the surface density as a function of fluence deduced from interference fringes. The fit was performed for fluences greater than 6 × 1011 ions cm−2 since the very beginning of the irradiation corresponds to ice compaction (e.g. Palumbo 2006; Dartois et al. 2013).
In the text

	Fig. 5. QMS signals as a function of fluence. Top panel: Background corrected signals. Bottom panel: Molecular fractions in the gas (solid lines, left y-axis) derived from the fragmentation patterns and ionisation cross sections, and sputtering yield of each species (dashed lines, right y-axis). See text for details.
In the text

	Fig. 6. Mass fragmentation pattern monitored during species injection (for ice sample deposition) for H2CO (top panel) and CO (unpublished data) (bottom panel).
In the text

Fig. 7. Schematic timeline of events during the first impact of a swift heavy ion (SHI) on an ice sample, within the thermal spike model. Chemistry and sputtering occur almost simultaneously within 10−9 s of the ion’s passage (Lang et al. 2020). The length ds corresponds to the sputtering depth probed by an individual ion incident, which is inferior to d0, the initial ice film thickness; rr corresponds to the radiolysis destruction radius, which is at most equal to the sputtering cylinder radius rs (Dartois et al. 2021). Considering the timescales involved, infrared spectroscopy cannot observe the radiolysis effects within the sputtered volume during the passage of an ion.

In the text

	Fig. 8. Comparison between the observed integrated intensities (blue dots) and the integrated intensities derived from the non-LTE model (red diamonds), for each of the observed lines, labelled by their rest frequencies. In the figure, ‘o’ stands for ortho and ‘p’ stands for para.
In the text

	Fig. 9. Intact sputtering rate of pure H2CO and effective sputtering rates of the pure ices CO, CO2, and H2O as a function of the ζ cosmic-ray ionisation rate. Effective sputtering yield pre-factors are taken from Dartois et al. (2023).
In the text

	Fig. 10. H2CO gas phase abundance with respect to H nuclei as a function of visual extinction. The observational value derived in this study for the prestellar core L1689B (blue hatched zone) is compared to five different models for H2CO sputtering: pure H2CO ice and H2CO ideally diluted in H2O, CO2, and CO matrices. See text for details.
In the text

	Fig. B.1. Observed H213CO spectra. The vertical dotted line shows the source velocity (v = 3.5 km s−1). The spin symmetry is indicated in the upper left corner of each plot, with ‘o’ standing for ortho and ‘p’ standing for para.
In the text

	Fig. B.2. Observed H2C18O spectra. The vertical dotted line shows the source velocity (v = 3.5 km s−1). The spin symmetry is indicated in the upper left corner of each plot, with ‘o’ standing for ortho and ‘p’ standing for para.
In the text

	Fig. C.1. Infrared transmittance spectra of a CO2 ice film evolution upon 74 MeV 86Kr18+ swift heavy ion irradiation. Solid lines: Experimental data. Dashed lines: Model spectra fitted to the data, without considering the absorption bands. The QMS signal ends around a fluence of 4.1012 ion cm−2, before the end of the irradiation as the instrument stopped prematurely.
In the text

	Fig. C.2. Evolution of the bulk mass per area as a function of fluence deduced from interference fringes for a CO2 ice film.
In the text

	Fig. C.3. QMS signals as a function of fluence. Top panel: Background corrected signals. Bottom panel: Molecular fractions in the gas (solid lines, left y-axis) derived from the fragmentation patterns and ionisation cross sections, and sputtering yield of each species (dashed lines, right y-axis). See text for details.
In the text

	Fig. C.4. Evolution of the background-corrected IR spectrum centered on the ν3 stretching mode of CO2 upon 74 MeV 86Kr18+ swift ion irradiation
In the text

	Fig. C.5. Column density derivative with respect to fluence. The solid and dashed lines are fits with or without including the sputtering depth Nd in Eq. C.3, respectively.
In the text