© The Authors 2026

1 Introduction

Molecules in the cold regions of the interstellar medium (ISM) are in their ground electronic and vibrational states and populate low-lying rotational states as temperatures in this environment range between 10 K and 100 K. Moreover, rotational spectroscopy has proven to be an important tool for detecting molecules in space (Herbst 2013; Weaver 2019). Analysis of the obtained molecular spectra can aid in understanding the evolution of cold cores of dense interstellar clouds into stellar and planetary systems. To accurately characterise the source of these spectra, a thorough and detailed understanding of the many elementary chemical-physical processes involving gaseous molecules and ions is necessary, both in their interactions with each other and with grain surfaces (Herbst 2013). Modelling these complex homogeneous and heterogeneous reactions at a microscopic level requires state-to-state resolved rate coefficients (Klessen & Glover 2016), particularly for rotational transitions, as the energy differences between rotational levels at low temperatures can play a fundamental role in determining whether a given reaction occurs. On the one hand, the advent of sophisticated experimental apparatuses allows scattering studies to be performed with rotationally controlled molecular beams (Toscano 2024); on the other hand, a complete theoretical evaluation of the transition probabilities between vibrational and rotational molecular states needs to achieve a proper description of the reaction dynamics. In recent years, we have developed a computational set-up useful for determining the scattering probabilities of molecules from surfaces resolved in vibrational and rotational states. This approach utilises accurate interaction potentials that explicitly consider long-range interactions between involved partners. These latter interactions play a crucial role since they drive, especially at low collision energies (E_coll), the formation of the precursor state with which each molecule, interacting with the surface, selectively evolves towards the final state of the process (see e.g. Rutigliano & Pirani 2021, 2020 and references therein). Molecular dynamics (MD) simulations are carried out using a semi-classical method that explicitly includes the coupling between the energy of impinging species and surface phonons (Billing 2000). Therefore, our chemical state-to-state MD simulations can be considered useful for investigating the occurring probability for elementary processes promoted by molecules impinging with low E_coll on surfaces frozen at low temperatures, thereby mimicking the conditions of the ISM. In particular, we can also obtain the adsorption (sticking) probability on the cold surface. This work focuses on the interaction of carbon monoxide (CO) with a graphite surface under conditions similar to those of the ISM. This choice stems from the fact that CO, present in appreciable concentrations in the ISM, participates in several elementary processes occurring in such an environment and its sticking coefficient on sub-micrometric cosmic grain analogues has been recently determined in the laboratory (Stadler et al. 2024). Furthermore, CO represents one of the reservoirs of carbon in space (Toscano 2024) and, being infrared active because of its heteronuclear character that is associated with a permanent electric dipole, it can be easily detected in galaxies, and hence it is exploited to estimate the masses of interstellar molecular clouds (Labiad et al. 2022). Precisely because of its peculiarities, CO is also often used in several environments, including the ISM, as a tracer of IR inactive diatomic molecules, mainly for the most abundant H₂ (Labiad et al. 2022). In order to obtain results that are of general interest, we also investigated, under the same conditions, the behaviour of N₂, a homonuclear (IR inactive) molecule that is the main component of air. We highlight the differences and similarities between CO and N₂ molecules, which are iso-electronic molecules with the same mass and very similar behaviours of their internal degrees of freedom.

2 Computational details

Molecular dynamics simulations have been carried out in the framework of a semi-classical collisional method, where the motion of the impinging molecule is treated classically by solving the corresponding Hamilton equations. The surface is treated quantum mechanically by solving the time-dependent Schrödinger equations of motion, under the harmonic oscillator approximation for the atoms in the surface perturbed by the force exercised by the impinging species. This method has been successfully used in the past (Rutigliano et al. 2001) to simulate, under low-energy flux, the recombination of atomic hydrogen on graphite, thus allowing the determination of values of probabilities for H₂ formation comparable to those obtained with quantum treatments, but with the added value of being able to resolve the roto-vibrational distributions of the formed molecules. The latter quantities are of considerable importance for studying the many processes occurring in the ISM. More recently, the method properly included, exploiting the improved Lennard-Jones function (Pirani et al. 2008), the role of long-range interactions that promote the formation of the physisorption potential well, whose features are crucial to control the balance between adsorption and scattering at the gas-surface interface (Rutigliano & Pirani 2016). The graphite surface model and, consequently, the eigenvalues and eigenvectors for the phonon motion assumed in the treatment are those of Rutigliano et al. (2001). In the present MD simulations, we use the same potential energy surfaces as Rutigliano & Pirani (2021) and Rutigliano & Pirani (2020). We set the surface temperature (T_S) at 10 K, and most of the attention is focused on CO molecules impinging on the surface in their ground roto-vibrational state (j_i, v_i) with E_coll ranging between 0.001 eV and 0.1 eV. The behaviour predicted for N₂ is used for a useful comparison. Selecting an approach normal to the surface plane, we propagate and examine 30 000 trajectories for each investigated E_coll value. The initial coordinates of the two atoms in the molecule are chosen randomly in an aiming area on the surface that is large enough, and at the same time able to prevent edge effects during trajectory propagation. The integration step is 0.25 fs and the accuracy required in the integration procedure is 10⁻⁸. The criteria adopted in the analysis of trajectories, leading to the assignment of a given trajectory to one of the possible reaction channels, are the same as those used in Rutigliano & Pirani (2021). Accordingly, molecular scattering occurs if, after the interaction with the graphite surface, the distance between the C and O atoms (r_C-O) in the molecule is smaller than the distance considered critical for CO molecule dissociation (r_diss ≈ 3.5 Å) and at the same time the distance between the molecule centre of mass (CM) and the surface is larger than 8.0 Å. On the contrary, if, after the interaction with the surface, the distance between the molecule CM and the surface is smaller than 5.0 Å and r_C-O < r_diss, the molecule is considered trapped in the physisorption well, while the molecule is dissociated if r_C-O > r _diss. The C and O atoms produced in the dissociation can be adsorbed or re-emitted in gas-phase as separate atoms; alternatively, one atom can remain adsorbed on the surface and the other can be re-emitted in gas-phase in the adsorption-desorption process (Rutigliano & Pirani 2021). The maximum propagation time for the trajectories is very long and fixed at 30 000 fs. From the trajectory analysis, we obtain the probability for the complete set of elementary processes occurring at the gas-surface interface and the final roto-vibrational states (j_f, v_f) in the case of scattered molecules. The probability of each elementary process is obtained by dividing the number of trajectories that produce that process by the total number of trajectories. The adopted method treats the molecular translational motion classically, whereas the rotational and vibrational populated levels of each scattered CO molecule, assumed to behave as a Morse oscillator (Huber & Herzberg 1979), are analysed in terms of the action-angle variables using the semi-classical quantisation rule (Muckerman 1971). According to this rule, at each trajectory time step, the internal energy of the interacting molecule is properly partitioned. In particular, the rotational energy is determined through relations, including rotational angular momentum and related spectroscopic constants, from which the corresponding rotational level is derived. The vibrational contribution is obtained as the difference from the total internal energy. In general, by averaging the occurring probability for a given process over a Maxwell-Boltzmann energy distribution evaluated at the temperature of interest, the corresponding effective coefficient is obtained. Accordingly, by averaging over the same energy distribution the state-to-state probability for CO scattering in a given rotational or vibrational (j/v) state, we can determine the corresponding state-selected reflection coefficient, Rj (R_v), for that state. In addition, energy accommodation coefficients, α_K, for the internal degrees of freedom of the molecule are also useful quantities (Klessen & Glover 2016). These coefficients, according to Goodman & Wachman (1976), are defined as α_K = (T_K - T′_K)/(T_K - T_S), where label K can be T (for translation), R (for rotation), and V (for vibration), and T_K and T′K are the temperatures corresponding to the given molecular degree of freedom before and after the interaction with the surface (taken at T_S), respectively. In this investigation we extended the analysis, carried out with high detail for CO impinging on the surface in its ground roto-vibrational state (j_i=0, v_i=0), also to the first excited rotational level (j_i =1, v_i=0), since the population sum of these two low-lying states falls in the range 70-100% for temperatures comparable to or lower than 10 K (typical of the ISM).The much higher state (j_i=5, v_i=0) is also considered for a useful comparison of the observables.

3 Results and discussion

The probability of the two occurring surface processes, namely molecular adsorption (sticking) and scattering for CO (0,0) and N₂ (0,0) colliding under thermal and sub-thermal conditions with the graphite surface is reported in Fig. 1 as a function of E_coll. This figure shows that CO remains adsorbed with high probability on the surface for collision energy lower than 0.01 eV. As E_coll increases, the probability of adsorption decreases, while that of scattering begins to rise.

An analysis of the trajectories reveals that CO interacts strongly with the surface and is preferentially captured when its C-end faces the surface. In contrast, when the O-end approaches the surface, molecules tend to be pushed away. This behaviour, highlighted in Rutigliano & Pirani (2021), is attributed to the torque effect. This effect, absent in the case of N₂ (Rutigliano & Pirani 2020), arises from the different strengths with which C and O interact with the surface, emphasised by the feature that the two atoms of different masses are placed asymmetrically with respect to the centre of mass of the molecule. Thus, anisotropic effects are enhanced for low collision energies and low surface temperatures. In Fig. 2, the probabilities of adsorption and scattering when CO impinges on the surface in the ground vibrational state and in two excited rotational states, j_i =1 and j_i=5, are reported as a function of E_coll in comparison with the same quantities obtained for the roto-vibrational ground state ( j=0, v_i=0).

The sticking coefficient determined in our simulation conditions results in 0.97, 0.96, and 0.58 for CO molecule impinging in j_i = 0, j_i = 1, and j_i = 5, respectively. These values, higher than those experimentally measured in Stadler et al. (2024) for low temperatures on grains, are in agreement with the values reported for an olivine crystal and reflect the statement that when increasing the size of the grain, the sticking coefficient increases to the limit (very close to 1) for flat and cold surfaces, as also emphasised for the interaction of H₂O and CO₂ under similar conditions (Laffon et al. 2021). The latter behaviour is also closely related to the presence of different sites on a larger surface area compared to that of sub-micrometre-sized grains. Moreover, it is important to consider that, in the experimental measurements, the molecules are not in well-defined rotational states. In Stadler et al. (2024), a comparison between CO and N₂ is also made. Accordingly, we also consider the interaction of N₂ in the roto-vibrational ground state, assuming the same simulation conditions as for CO, and obtained results for the adsorption and scattering probability are reported in Fig. 1 again for a useful comparison. The adsorption of N₂ is less probable than that for CO, and its probability decreases approximately linearly with the collision energy. The sticking coefficient calculated for N₂ is 0.71, lower than that for CO, in agreement with the previous result (Stadler et al. 2024). However, in Stadler et al. (2024) the sticking coefficients measured on sub-micrometre grains are 0.05±0.02 and 0.07±0.02 for CO and N₂, respectively. These values are nearly identical, within the error, suggesting that on sub-micrometre substrates and under low-temperature conditions the two molecules behave similarly and that the interaction dynamics is governed only by the molecular mass. Therefore, the presence of a larger and flat surface at the same temperature, on the other hand, is essential to highlight the different behaviours of the two molecules at low collision energies and low surface temperatures. Consequently, the picture of interaction dynamics on larger grains (i.e. exceeding sub-micrometre sizes) can be quite different.

Molecular dynamics simulations indicate that impinging CO is captured in a region near the surface and, as previously emphasised, becomes trapped in a roto-translational motion, predominantly governed by the end facing the surface, as previously said. Additionally, the amplitudes of the oscillations occurring on the surface are smaller compared to those observed for the N₂ molecule. A typical trajectory that ends with CO adsorption and reflects this behaviour is shown in Fig. 3, where the time evolution is reported along the three spatial directions for a typical adsorption trajectory, thus showing the torque effect.

In Rutigliano & Pirani (2021), we reporte values of 0.771 and 5.04 ×10⁻³ for the CO adsorption probability at E_coll=0.01 eV and E_coll=0.1 eV, respectively, when the surface is taken at T_S=100K. Comparing the corresponding values for the coldest surface shown in Fig. 1, it can be inferred that as the surface temperature increases, the adsorption decreases slightly. The observed decrease may be related to an increase in thermal desorption, which is primarily noticeable for CO, which interacts more strongly, particularly with the C-end, with the surface. In contrast, for N₂, the influence of surface temperature on the calculated probabilities for surface processes falls well within the expected error margin of 10% for the values derived from the adopted simulation methodology.

The stereodynamic effects in the energy transfers observed during CO adsorption are also evident in the rotational distributions of scattered CO molecules, while these effects are missing in the corresponding distributions of N₂ molecules scattered under the same conditions. In Fig. 4, the final rotational distributions achieved by the two molecules, impinging on the surface in their roto-vibrational ground-state (v_f=0, j_i =0) at E_coll=0.001 eV and E_coll=0.01 eV and T_S = 10 K, are reported. This figure shows that the major anisotropic character of the interaction with the surface of heteronuclear CO, whose role mainly emerges at the lowest probed Ecoll, produces wider distributions, thus populating low-lying to medium-lying final rotational levels (j_f) between j_f=0 and j_f=7; instead by increasing E_coll, the distribution peaks at j_f=1. On the contrary, for N₂, the determined distributions are very similar for the two collision energies. This behaviour is ascribable to the circumstance that the two ends of the molecule undergo an interaction of the same strength with the surface.

The distribution obtained for CO molecules scattered at the lowest Ecoll could be thought of as the result of a kind of resonance established between the collision energy with which the molecule approaches the surface, corresponding to about 11.6 K in terms of temperature scale, and the temperature of the surface, influencing the motion of the latter’s phonons. This resonance pumps energy into the CO molecule, populating higher-lying rotational levels, whilst the molecules leave the surface with higher translational energy than when they arrived, as assessed by determining the corresponding energy accommodation coefficient (see below). What has been described does not occur in the case of N₂, suggesting that this phenomenon is a specific peculiarity of the CO molecule. This effect can be attributed entirely to the greater interaction of the C-end of the molecule with the surface, and to the torque effect (Rutigliano & Pirani 2021). As a result, based on the definition of the energy accommodation coefficient, the coefficients for translation and rotation are α_T = −6 and α_R = 9.4, respectively. On the contrary, for E_coll=0.01eV, the greater mismatch between the surface temperature and the collision energy does not produce resonance, and only the level j_f =1 populates with a probability of ~43% (see Fig. 4). The remaining part of the available energy is funnelled into the translational motion of the molecule, which consequently leaves the surface with lower translational energy than when it arrived. Accordingly, the energy accommodation coefficient for the translational motion is only 0.8. For the initial rotational state j_i=0, the determined state-to-state reflection coefficient, Rj, promoting the formation of final states, j_f, is reported in Table 1. From this Table we can deduce that the highest coefficient corresponds to j_f=3.

The same analysis was also carried out for CO impinging on the surface with the other low-lying roto-vibrational state (j_i = 1, v_i=0), significantly populated in the ISM. It has been found that the small increase in molecular internal energy (of the order of a Kelvin) does not substantially modify the resonance features that promote rotational excitation of desorbed molecules (see Fig. 5). We also investigated the behaviour of CO with initial state (j_i=5, v_i=0), whose internal energy here amounts to about 10² K. We find that under such conditions the related collision events promote CO molecules scattered in rotational levels having a different behaviour than that observed for the lowest considered j_i=0 value. A comparison of Figs. 4 and 5 shows, for j_i=1, that the distribution for E_coll=0.001eV preserves the main features observed for j_i = 0, even if slightly attenuated; instead, for j_i = 5 the distribution behaves similarly for both E_coll = 0.001 eV and E_coll=0.01 eV. Moreover, at low E_coll, additional information on the selectivity of the scattering process is provided in Table 1, where the state-resolved reflection coefficients, R_j, for the most populated levels, for the CO in the different j_i values considered here, are reported for comparison.

These findings have important implications since CO molecules, rotationally excited by the resonance with the surface phonons, remain in the gas phase conserving their excitation. They then surely affect the IR response of a detector that adopts CO as a tracer of other IR inactive species in cold environments such as the ISM.

Our surface model could be considered as a precursor of polycyclic aromatic hydrocarbons (PAHs) or of carbon grains. The radii for PAHs range between 5 Å and 10 Å, while for standard grains the radii range between 300 Å and 3.0 × 10⁴ Å (Herbst et al. 2005). Assuming a round grain, we could estimate the cross-section for the scattering of CO (0,0) as a function of E_coll. In the limit of the low extremity value for the radius for PAHs and grains, 5 Å and 300 Å, respectively, we obtain, at E_coll=0.001eV, a cross-section of 2.51 Ų for PAHs and 9.05 × 10³ Ų for the grains. From the cross-section by integrating on a Maxwell-Boltzmann distribution for the energy, we obtained the resolved state-to-state rate constant for the most populated states, as these quantities are particularly useful in microscopic models (Klessen & Glover 2016). As an example, the state-to-state rate constants for the three initial values of the rotational states here considered are reported in Table 2.

It is worth noting once again the peculiar behaviour of CO in the two lowest rotational states, triggered by a more favourable condition for the energy exchange between molecule and substrate.

Fig. 1 Probability for molecular adsorption and scattering for CO (red symbols) and N2 (blue symbols) as a function of collision energies and for TS=10 K. The grey lines are included for clarity.

Fig. 2 Probability for molecular adsorption and scattering for CO (0,1), CO (0,5), in comparison with the same quantities for CO (0,0), as a function of collision energies and for TS=10K. The grey lines are included for clarity.

Fig. 3 Motion in X, Y, Z for C (blue line) and O (red line) for a typical adsorption trajectory for CO (0,0) impinging the surface with Ecoll=0.001 eV. TS=10 K.

Fig. 4 Normalised final rotational distributions for CO (0,0) (in red) and N2 (0,0) (in blue) impinging the surface with Ecoll=0.001 eV (upper panels) and Ecoll=0.01 eV (lower panels) scattered from a graphite surface taken at 10 K.

Fig. 5 Normalised final rotational distributions for CO (0,1) (in blue) and CO (0,5) (in red) impinging the surface with Ecoll=0.001 eV (upper panels) and Ecoll=0.01 eV (lower panels) scattered from a graphite surface taken at 10 K.

4 Conclusions

As CO collides with H₂ in cold environments, it experiences rotational excitation in low-lying levels followed by radiative decay. These energy transfers are some of the most commonly used methods to trace molecular hydrogen, which is invisible to radio telescopes in the dense ISM. In particular, state-to-state rate coefficients for rotational energy transfer of CO (j_i=0, 1, 4 → j_f =0–3) in collisions with H2, measured and calculated at low temperature (T= 20, 10, and 5.5 K), are quite high (for full details, see Labiad et al. 2022). From our results it emerges that in the low collision energy regime, in a low-temperature environment such as the ISM, rotationally excited CO molecules can also be generated by the gas-surface interaction, such as that controlling adsorption–desorption on large graphitic grains. In particular, it has been found that a sort of resonance between the low energy of molecules, impinging in the lowest roto-vibrational state, and that of the phonons of the cold surface promotes the formation with a significant probability of excited rotational states between j_f =1 and j_f=6. Furthermore, if these rotationally excited molecules collided again with the grain surface, they would most likely remain excited, as evidenced by the calculations done with j_i =1 and j_i=5. Therefore, from our results, it emerges that the stimulated CO rotational radiative decay cannot be unequivocally ascribed to the collisions of CO with H₂, as CO rotationally excited in the same levels can also be directly produced by adsorption-desorption on the cold grain’s surface, although state-to-state rate coefficients are some orders of magnitude smaller. We note that the difference in the rotational excitation selectivity of CO and N₂ (which are isoelectronic molecules with the same mass), promoted by their collisions at low collision energy with a cold graphite surface and of interest for the stereodynamics of the basic elementary processes occurring at the gas-surface interface, tends to disappear at high temperatures, and this topic is of interest for fundamental elementary processes occurring under hyper-thermal conditions in the gas phase. In particular, the characterisation of thermochemical non-equilibrium effects in gaseous mixtures heated to some thousands of Kelvin is one of the research frontiers in hypersonic flight of vehicles or spacecraft re-entering the Earth’s atmosphere (He et al. 2024). Current attention is focusing on the possibility of controlling the vibrational excitation behaviours of shock-heated high-temperature air, mostly formed by N₂ (about 80%) and O₂ (about 20%), which are both IR inactive, by using the CO roto-vibrational thermometry. A recent combination of experimental results with theoretical simulations highlights that CO cannot be used as a trace gas to detect the thermal nonequilibrium effects of high-temperature O₂ (He et al. 2025) since it exhibits a different vibrational temperature with respect to CO. On the other hand, the vibrational temperature of the CO trace in the air must exhibit time histories consistent with those of N₂, its main component, because of the efficient vibrationvibration energy transfer between CO and N₂, and this study is still in progress. In conclusion, it is emerging that a deepened and extended characterisation of the CO behaviour under various conditions, covering sub-thermal to hyper-thermal regimes, opens an interesting perspective of applications in several areas of fundamental and applied research.

References

All Tables

All Figures

	Fig. 1 Probability for molecular adsorption and scattering for CO (red symbols) and N2 (blue symbols) as a function of collision energies and for TS=10 K. The grey lines are included for clarity.
In the text

	Fig. 2 Probability for molecular adsorption and scattering for CO (0,1), CO (0,5), in comparison with the same quantities for CO (0,0), as a function of collision energies and for TS=10K. The grey lines are included for clarity.
In the text

	Fig. 3 Motion in X, Y, Z for C (blue line) and O (red line) for a typical adsorption trajectory for CO (0,0) impinging the surface with Ecoll=0.001 eV. TS=10 K.
In the text

	Fig. 4 Normalised final rotational distributions for CO (0,0) (in red) and N2 (0,0) (in blue) impinging the surface with Ecoll=0.001 eV (upper panels) and Ecoll=0.01 eV (lower panels) scattered from a graphite surface taken at 10 K.
In the text

	Fig. 5 Normalised final rotational distributions for CO (0,1) (in blue) and CO (0,5) (in red) impinging the surface with Ecoll=0.001 eV (upper panels) and Ecoll=0.01 eV (lower panels) scattered from a graphite surface taken at 10 K.
In the text