© The Authors 2026

1 Introduction

Comets are thought to be remnants of planet formation and preserve valuable clues about the origin of our Solar System (Mumma & Charnley 2011). These objects may have played a role in delivering organic molecules that acted as building blocks of life on Earth. The nucleus of a comet is composed of rock surrounded by ices, consisting primarily of H₂O but also of CO and CO₂. As the comet approaches the Sun, these ices sublimate due to increasing solar radiation, forming an atmosphere around the nucleus known as the cometary coma. For comets at short heliocentric distances (r_h ~ 1 a.u.), the surface temperature driven by solar radiation exceeds the H₂O sublimation temperature, making water the dominant species in the coma.

In the coma, the density goes from 10¹² cm⁻³ close to the surface of the nucleus to 1 cm⁻³ in the outer coma (Despois et al. 2006). At intermediate nucleocentric distances, the density is neither high enough for local thermodynamic equilibrium (LTE) conditions to prevail nor low enough for the fluorescence regime to be achieved. This can lead to strong non-LTE signatures in the molecular lines (Bodewits et al. 2022). Hence, the competition between radiative and collisional (de-)excitation processes must be accounted for, and thus collisional rate coefficients are required. As most cometary observations are performed on comets at short heliocentric distances (r_h ~ 1 a.u.), where the molecular outgassing from the comet is the highest, the collisional rate coefficients with H₂O as a collisional partner are therefore primarily required.

The determination of collisional rate coefficients is recognized as a challenging task, and it is mostly carried out using numerical calculations involving rare gases and light molecules such as H₂ for interstellar applications. Since water has a relatively high mass compared to these colliders and possesses a large dipole moment, the number of channels required to achieve convergence in the calculations makes the use of quantum approaches unfeasible. Alternative methods are necessary. Currently, among over 50 molecules detected in comets through remote observations (Biver et al. 2022), collisional data with the H₂O projectile exist for four molecules only. The most recent sets of collisional rate coefficients are CO–H₂O by Loreau et al. (2018a); Faure et al. (2020), HF–H₂O by Loreau et al. (2020, 2022), H₂O–H₂O data by Mandal & Babikov (2023b,a), and HCN–H₂O by Żółtowski et al. (2025).

In these studies, the scattering calculations were performed using two scattering approaches: the mixed quantum classical trajectory [MQCT, Semenov & Babikov (2013); Mandal et al. (2024)] for the H₂O–H₂O data, and the statistical adiabatic channel model [SACM, Loreau et al. (2018b)] approaches for the other three systems. Both approaches are complementary: MQCT relies on classical mechanics and is therefore inherently more accurate at high temperatures. It has demonstrated an acceptable accuracy down to 100 K according to the work of Joy et al. (2024). SACM relies on statistical theory and is thus inherently more accurate at low temperatures. It has been shown to provide accurate results up to 100 K by Loreau et al. (2018a) in their study of CO–H₂O rate coefficients.

The SACM approach employed in this work was initially proposed by Quack & Troe (1974, 1975) to treat inelastic and reactive collisions. It assumes that the complex formed by the two colliders lives long enough for the total energy of the collision to be statistically redistributed. The method operates under the straightforward assumption that all energetically accessible channels¹ have equal probabilities, while all others have zero probability. Loreau et al. (2018b) refined this approach by accurately computing the adiabatic states corresponding to each channel, rather than relying on simple analytical formulas.

CS is one of the most abundant sulfur-bearing species detected in comets (Biver et al. 2022), and it has been clearly identified as a daughter molecule – meaning that it is produced by the photolysis of a species directly outgassing from the nucleus (Roth et al. 2021b). However, its parent molecule remains unidentified. CS₂ was initially proposed as a possible parent molecule, but it was later found to have a lifetime that is too short and an abundance that is too low for it to be at the origin of CS formation (Calmonte et al. 2016; Roth et al. 2021b; Biver et al. 2022). Recently, Biver et al. (2024) reported that the photodissociation rate of H₂CS closely matches the rate previously estimated for the parent molecule of CS, supporting the hypothesis that H₂CS could be a plausible parent of CS in cometary comae. However, H₂CS is also too scarce to be a likely parent. Therefore, no suitable candidate for the parent molecule of CS currently exists.

In this work, we determined the collisional rate coefficients characterizing the inelastic collisions $\mathrm{C}\mathrm{S}\left(j_{\mathrm{C}\mathrm{S}}\right)+{\mathrm{H}}_2\mathrm{O}\left(j_{k_a{k}_c}\right)\to \mathrm{C}\mathrm{S}\left(j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}\right)+{\mathrm{H}}_2\mathrm{O}\left(j_{k_a^{\mathrm{\prime }}{k}_c^{\mathrm{\prime }}}^{\mathrm{\prime }}\right)$ from 5 to 100 K for cometary applications. For that purpose, the CS–H₂O potential energy surface (PES) characterizing the electronic interaction between the colliders was computed using the SAPT-DFT. The scattering calculations were performed based on the computed PES using the SACM approach. The rate coefficients were computed by including all transitions between CS levels up to j_CS ≤ 18 (279.44 cm⁻¹), and para-H₂O (hereafter, noted p-H₂O) levels up to the j_kakc ≤ 2₂₀ (136.56 cm⁻¹), and ortho-H₂O (hereafter, noted o-H₂O) levels up to j_kakc ≤ 3₀₃ (136.91 cm⁻¹). The new CS–H₂O of rate coefficients have been employed in non-LTE excitation models, alongside CS–H₂ rate coefficients and the rate coefficients derived by Biver et al. (1999). The latter are derived based on a global cross section. The objective is to assess the impact of the different datasets on the excitation conditions of CS in comets. This paper also includes a discussion of the expected accuracy of the CS–H₂O data and proposes guidelines for which approximations can be used to accurately model the excitation of molecules in comets when collisional rate coefficients with H₂O are unavailable. Finally, conclusions are presented.

2 Numerical methods

2.1 Ab initio calculations

We developed a new ab initio PES for the scattering calculations between the CS and H₂O colliders. We considered CS and H₂O molecules as rigid with their geometries corresponding to the ground electronic and rovibrational states. The CS interatomic distance was fixed at 2.9179a₀ and r(OH) bond distances in the H₂O molecule were fixed at 1.8368a₀ with the HOH angle of 104.7 degrees. This leads to the five-dimensional PES; namely, the interaction is completely described by one radial (R) and four angular coordinates.

We followed the recent work of Dagdigian (2020) to set up the origin of the coordinate system in the center of mass (c.o.m.) of the H₂O monomer. With this choice, the V(R, θ₁, ϕ₁, θ₂, ϕ₂) PES describes the position of the c.o.m. of the CS molecule around the H₂O monomer with a Jacobi vector, R, which has (R, θ₁, ϕ₁) spherical coordinates. (R, θ₁, ϕ₁) has its origin on the H₂O c.o.m., with R being the intermolecular distance between the c.o.m.s of H₂O and CS and the pair of angles (θ₁, ϕ₁) describing the angular position of the c.o.m. of the CS molecule. The remaining spherical angles (θ₂, ϕ₂) describe the relative orientation of the CS rotor with respect to an H₂O fixed system of coordinates, which has a z axis coincidental with the C₂ symmetry axis of the water molecule and in which the water molecule lies on the xz plane.

To determine the PES, we used the SAPT(DFT) (Misquitta & Szalewicz 2005; Misquitta et al. 2005) with the augmented correlation-consistent quadruple-zeta basis set (aug-cc-pVQZ). The PES was generated using the AUTOPES package (Metz et al. 2016). In this approach, the short-range part of the interaction is obtained by partitioning the system into two monomers (H₂O and CS in our case) and evaluating the electrostatic, exchange, induction, and dispersion energies together with their exchange counterparts (exchange–induction and exchange–dispersion). In addition, the $\delta E_{\text{int,resp}}^{\mathrm{H}\mathrm{F}}$ correction (Patkowski et al. 2006) was included in the final SAPT-DFT PES to account for higher-order induction and exchange–induction effects. The long-range part was constructed independently from monomer electron densities and frequency-dependent density susceptibility functions at the coupled Kohn–Sham level, leading to a c.o.m.–c.o.m. asymptotic multipole expansion of the interaction energy in inverse powers of R. The two parts are then smoothly merged into a global PES, which is built by freezing the damped asymptotic part and fitting only the short-range term, which produces a smooth potential (Metz et al. 2016). The resulting surface reproduces the interaction energies with root-mean-square (RMS) errors of 2.5 cm⁻¹ for energies below −600 cm⁻¹, 18 cm⁻¹ along the repulsive wall up to 3500 cm⁻¹, 4.5 cm⁻¹ for all bound configurations (energies below 0), and 57 cm⁻¹ for highly repulsive configurations near 35 000 cm⁻¹. Those errors are not expected to significantly impact the scattering calculations, considering that they are marginal with respect to both the well depth of the system and the accuracy of the SACM method.

Additionally, we performed benchmark calculations employing the “gold standard” of electronic structure methods, the coupled cluster method with single, double, and non-iterative triple excitations with explicit correlation [CCSD(T)-F12], using the augmented correlation-consistent triple-zeta basis set (aug-cc-pVTZ) supplemented with mid-bond functions for selected cuts of the PES. A good agreement was found between the SAPT(DFT) and CCSD(T)-F12 PESs. In Fig. 1 we show the comparison of the SAPT-DFT fit to selected cuts from our calculations with the gold standard CCSD(T)-F12 method. We note that the CCSD(T)-F12 calculations are more computationally demanding than SAPT(DFT) calculations, especially to scan five angular dimensions. For the relative orientation of CS and H₂O near the global minimum (geometry labeled “C” in Fig. 1), the CCSD(T)-F12 PES gives the well depth of 1269.6 cm⁻¹ at R = 8.07a₀, which is only about 50 cm⁻¹ deeper (4% relative difference) than the SAPT-DFT minimum. This difference in the global minimum well depth is mainly due to the inclusion of the perturbative triple excitations in the CCSD(T)-F12 method. We found by computing interaction energies without mid-bond functions that the bond functions contribute about 3.4 cm⁻¹ to the global minimum for the arrangement C shown in Fig. 1. For arrangements A and B, the mid-bond functions contribute on a similar order, about 5 cm⁻¹, which is about 1–2% of the interaction energy at minima of these arrangements. The CCSD(T)-F12 method with triple-zeta basis and mid-bond functions provides in practice interaction energies close to the complete basis set limit.

The next step after completion of the electronic structure calculations of the PES is to represent the potential as an expansion in the tensor spherical functions used in HIBRIDON program. We followed the work of Dagdigian (2020) for this step and we refer the reader to this reference for details of the angular expansion. The only difference is that Dagdigian (2020) considered H₂ as a collider, which limits the rotational quantum number, l₂, of the spherical rotor function in the potential expansion to being evenvalued, which is not the case for our CS molecule, where both spherical rotor functions with odd and even-valued l₂ are necessary. We expanded the PES to include expansion coefficients up to l₁ = 10, l₂ = 6, and l = 6. The expansion of the SAPT-DFT PES in those terms results in the following RMS errors: at the distance R = 5a₀, the RMS is 128 cm⁻¹ for interaction energies ranging from 10³ to 10⁴ cm⁻¹; while at R = 8.1a₀ (near the R_e of the PES), the RMS is 2.7 cm⁻¹ for energies ranging from - 1100 cm⁻¹ to 900 cm⁻¹. Finally, for the intermediate (R = 10a₀, and energy range from −500 cm⁻¹ to 300 cm⁻¹) and longer range (R = 20a₀, and an energy range from −30 cm⁻¹ to 30 cm⁻¹), the RMSs are 0.33 cm⁻¹ and 0.001 cm⁻¹, respectively. We conclude that this expansion provides a good description of the PES for our purpose without the significant increase in computational time that larger expansion sets may pose.

Fig. 1 Comparison of the CS–H2O SAPT(DFT) and aVQZ PES fit (solid black line) with the CCSD(T)-F12 and aVTZ+BF calculations performed for selected radial cuts. The θ1 and ϕ1 are spherical angles of the R vector, while θ2 and ϕ2 are spherical angles describing the orientation of the CS molecule.

2.2 Scattering calculations

To obtain all of the information about the collision, the scattering matrices, denoted as $S_{\gamma ,{\gamma }^{\mathrm{\prime }}}^J(E)({\gamma }^{\mathrm{\prime }\mathrm{\prime }}\equiv j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }\mathrm{\prime }}{j}^{\mathrm{\prime }\mathrm{\prime }}{k}_a^{\mathrm{\prime }\mathrm{\prime }}{k}_c^{\mathrm{\prime }\mathrm{\prime }}{j}_{12}^{\mathrm{\prime }\mathrm{\prime }}{l}^{\mathrm{\prime }\mathrm{\prime }})$, were calculated for each total angular momentum J. The SACM approach starts with the same equations then the close-coupling (hereafter, CC) approach, expressed as follows (Alexander et al. 2023): $$[\mathbf{I}\frac{d^2}{d{R}^2}+\mathbf{W}(R)]\mathbf{F}(R)=0,$$(1)

with the identity matrix, I, the radial part of the wavefunction, F(R), and the W(R) matrix defined as $$\mathbf{W}(R)={\mathbf{k}}^2-{\mathbf{l}}^2-\frac{2\mu }{{\hslash }^2}\mathbf{V}(R),$$(2)

with μ the reduced mass of the collisional system. V(R) is the matrix of the coupling potential, and k and l are diagonal matrix of the wavevector and of the relative angular momentum of the collisional partners, defined, respectively, as $$k_{ii}^2=\frac{2\mu }{{\hslash }^2}(E-{ϵ}_i),$$(3) $$l_{ii}^2=\frac{{\hslash }^2}{2\mu R^2}{l}_i(l_i+1),$$(4)

with ℏ the Planck constant divided by 2π, ϵ_i the energy of the channel, i, and l_i the relative orbital angular momentum of the i^th channel. The diagonalization of the V(R) + l² matrix yields the diagonal matrix of adiabatic wavevectors, which define the locally adiabatic states.

To study the inelastic CS–p-H₂O and CS–o-H₂O collisions using the SACM approach, we computed the adiabatic states of the system by diagonalization of the potential matrix, V(R), including the relative angular momentum matrix, l², following the work of Loreau et al. (2018a). A probability was then assigned to the channels – based on whether they are open or closed – to construct the $S_{\gamma ,{\gamma }^{\mathrm{\prime }}}^J(E)$ matrices. For a channel to be open, it must be energetically accessible at the total energy of the system, E. It means that its adiabatic state must present a minimum, and that its energy between the minima and its asymptotic value must always be lower than E. A closed channel is one that does not meet one of these requirements. The SACM approach works under the straightforward assumption that all open channels has the same probability, $\frac{1}{N}$, with N the total number of open channels, while any closed channel has zero probability. Additional details of the method can be found in Loreau et al. (2018b,a), Żółtowski et al. (2025), Godard Palluet et al. (2025), and Tonolo et al. (2025).

From these $S_{\gamma ,{\gamma }^{\mathrm{\prime }}}^J(E)$ matrices, the inelastic cross sections were finally computed as $${\sigma }_{\gamma \to {\gamma }^{\mathrm{\prime }}}(E)=\frac{\pi }{g_1{g}_2{k}_{\gamma ,{\gamma }^{\mathrm{\prime }}}^2}\sum _J(2J+1)\sum _{l{l}^{\mathrm{\prime }}}{|S_{\gamma ,{\gamma }^{\mathrm{\prime }}}^J(E)|}^2,$$(5)

with g_i the degeneracy of the rotational state of monomer i, and $k_{\gamma ,{\gamma }^{\mathrm{\prime }}}^2=\frac{2\mu }{{\hslash }^2}(E-E_{\gamma })$ the wavenumber.

The state-to-state rate coefficients at a given kinetic temperature, T_k, were obtained by calculating the thermal average of the cross sections over kinetic energies $$k_{\gamma \to {\gamma }^{\mathrm{\prime }}}\left(T_k\right)=\sqrt{\frac{8}{\pi \mu k_B^3{T}_k^3}}\int {\sigma }_{\gamma \to {\gamma }^{\mathrm{\prime }}}\left(E_k\right)E_k\exp \left(\frac{-E_k}{k_B{T}_k}\right)d{E}_k,$$(6)

where E_k = E − E_γ is the kinetic energy of the transition.

With this approach, the computational cost is drastically reduced for three reasons:

1. The determination of the adiabatic states does not require a R grid as thin or as large as for the propagation of F(R), which is the matrix defining the expansion coefficients. Indeed, only the region between the minimum and the centrifugal barrier region must be accurately described; the long-range interaction will not impact the adiabatic counting.

2. The convergence of the cross sections does not require the inclusion of asymptotically closed channels in the calculation of V(R) and l² matrix elements, contrary to quantum approaches such as CC or the coupled-state approximation. Thus, a much smaller rotational basis can be used, significantly reducing the computational cost.

3. The adiabatic states are independent of the total energy, E, so they are computed once per J. Only the adiabatic channel counting – the least expensive part of these calculations – needs to be performed at each energy.

The central processing unit (CPU) time required for both CC and SACM approaches is proportional to $N_{ch}^3$ (N_ch being the number of channels). Therefore, if N_ch is reduced by a factor of 5 with the SACM method, the CPU time will be reduced by a factor of 5³. An additional factor will also saved considering that the adiabatic states are valid at all energies.

In this work the goal was to compute rate coefficients characterizing the $\mathrm{C}\mathrm{S}\left(j_{\mathrm{C}\mathrm{S}}\right)+{\mathrm{H}}_2\mathrm{O}\left(j_{k_a{k}_c}\right)\to \mathrm{C}\mathrm{S}\left(j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}\right)+{\mathrm{H}}_2\mathrm{O}\left(j_{k_a^{\mathrm{\prime }}{k}_c^{\mathrm{\prime }}}^{\mathrm{\prime }}\right)$ inelastic collisions for temperatures from 5 to 100 K. The rate coefficients thus include transitions involving CS rotational levels of up to j_CS ≤ 18, and the first five levels of p– and o-H₂O, up to j_kakc ≤ 2₂₀ for p-H₂O, and up to j_kakc ≤ 3₀₃ for o-H₂O.

The rotational constants of water were taken from Herzberg (1966), with A = 27.877 cm⁻¹, B = 14.512 cm⁻¹, and C = 9.285 cm⁻¹ the rotational constants corresponding to the moments of inertia along the a, b, and c axes, respectively, with I_a ≠ I_b ≠ I_c. The rotational constant of CS was taken as B_CS = 0.817 cm⁻¹ according to Bustreel et al. (1979).

The adiabatic states were computed using the HIBRIDON software (Alexander et al. 2023). To reduce the computational cost, the number of channels was limited by using both the EMAX parameter and the rotational basis. The EMAX option will remove from the basis any channels with an energy greater than EMAX². As the convergence of the cross sections is weakly sensitive to asymptotically closed channels within the SACM approach, using this option will not impact the precision of the results, but will reduce the computational time. The rotational basis was finally restricted to j_CS ≤ 31 and j_kakc ≤ 3_kakc for CS–p-H₂O, and to j_CS ≤ 26 and j_kakc ≤ 3_kakc for CS–o-H₂O to reduce computational cost while maintaining the precision of the rate coefficients. The EMAX parameter was set at 817 cm⁻¹, the highest total energy considered in the scattering calculations. The convergence of the rotational basis and of the EMAX parameter for the adiabatic state calculations, along with additional information about the calculations, is presented in Appendix A.

The cross sections were then calculated following the SACM approach described above, and summed over an increasing number of partial waves, J, until a convergence of the cross sections at better than 5% was observed. As a result, partial waves with 0 ≤ J ≤ 153 were included for the calculation of CS–p-H₂O cross sections, and with 0 ≤ J ≤ 149 for the CS–o-H₂O cross sections. By averaging the cross sections over a thermal distribution of the kinetic energies, the state-to-state rate coefficients were obtained from Eq. (6).

2.3 Non-LTE modeling

To evaluate the influence of the new rate coefficients on the interpretation of cometary observations, we incorporated the CS–H₂O rate coefficients into a non-LTE excitation model to simulate the population of CS energy levels at various nucleocentric distances within the coma. The corresponding gas density was evaluated using the Haser (1957) model as follows: $$n_{{\mathrm{H}}_2\mathrm{O}}(r)=\frac{Q_{{\mathrm{H}}_2\mathrm{O}}}{4\pi r^2{v}_{exp}}\exp (-r\frac{{\beta }_{{\mathrm{H}}_2\mathrm{O}}}{v_{exp}}),$$(7)

with n_H2O the density of water vapor as a function of the nucleocentric distance, r, and β_H2O = 1.042 × 10⁻⁵ s⁻¹ the photodissociation rate of water taken from Zakharov et al. (2007). Q_H2O, the production rate of water, was set at 10²⁹ cm⁻³ s⁻¹. v_exp, the expansion velocity of the gas, was set at 1 km s⁻¹. The value of the production rates and the expansion velocity was settled in order to represent the “average” water-dominated comet at a heliocentric distance of r_h = 1 a.u. (Bodewits et al. 2022), following the methodology set out in Lupu et al. (2007) and Faure et al. (2020).

The population of CS energy levels at different nucleocentric distances was determined by solving the coupled radiative transfer equations and statistical equilibrium equations using the escape probability formalism as implemented in the RADEX code developed by van der Tak et al. (2007). The nucleocentric distances were considered between 10² and 10⁶ km. The kinetic temperature of the gas was set at 100 K.

In comets, CS represents between 0.03 and 0.2% of cometary volatiles, and CO represents on average about 10% (Biver et al. 2022). Therefore, the column density of CS, denoted as N_CS, was set at 1 × 10¹² cm⁻², so 1% of the CO column density reported by Lupu et al. (2007), who estimated a CO column density of approximately 1 × 10¹⁴ cm⁻² at a nucleocentric distance of r ≃ 10³ km in a water-dominated comet at an heliocentric of r_h ~ 1 a.u³.

The linewidth was set at δν = 2v_exp. The background radiation field was considered as the combination of the cosmic microwave background radiation field at 2.73 K and the solar radiation field. The latter was considered as a blackbody of 5770 K diluted with the corresponding factor W = Ω_S/4π, where Ω_S ≃ 6.79 × 10⁻⁵ sr is the solid angle of the Sun at a heliocentric distance, r_h, of 1 a.u. The Einstein coefficients were taken from the Cologne Database for Molecular Spectroscopy (Müller et al. 2001, 2005; Endres et al. 2016).

The RADEX code requires thermalized rate coefficients, which were computed from the state-to-state rate coefficients by assuming a thermal distribution of o/p-H₂O as $$k_{j_{\mathrm{C}\mathrm{S}}\to j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}}(T)=\sum _{j_{k_a{k}_c}}{n}_{j_{k_a{k}_c}}(T)\sum _{j_{k_a^{\mathrm{\prime }}{k}_c^{\mathrm{\prime }}}^{\mathrm{\prime }}}{k}_{\gamma \to {\gamma }^{\mathrm{\prime }}}(T),$$(8)

with n_jkakc (T) the population of the j_kakc level at the kinetic temperature, T: $$n_{j_{k_a{k}_c}}(T)=\frac{(2j+1)\exp \left(\frac{-E_{j_{k_a{k}_c}}}{k_{B}T}\right)}{\sum _{j_{k_a^{\mathrm{\prime }}{k}_c^{\mathrm{\prime }}}^{\mathrm{\prime }}}(2j^{\mathrm{\prime }}+1)\exp \left(\frac{-E_{j_{k_a^{\mathrm{\prime }}{k}_c^{\mathrm{\prime }}}^{\mathrm{\prime }}}}{k_{B}T}\right)}.$$(9)

The thermalized rate coefficients were used by assuming an o/p ratio of three, which is a typical value observed in comets at a kinetic temperature higher than 50 K (Biver et al. 2022).

This model is a very simplified one compared to the complex processes taking place in the coma: the chemistry of CS is neglected, as is the effect of the heliocentric distance on the abundance of CS, the impact of electrons, and the excitation of CS vibrational states by solar radiation. The objective of this model is to capture collisional effects, and we do not intend to predict any quantitative data from this model. For the accurate interpretation of observational molecular spectra of CS, more sophisticated approaches exist, such as the SUBLIME code described in Cordiner et al. (2023).

Fig. 2 Global minimum equilibrium geometry of the CS–H2O SAPT(DFT) and aVQZ PES. The Re between the c.o.m. of H2O and c.o.m. of CS is 8.11a0, with θ1 = 122°, θ2 = 297°, and ϕ1 = ϕ2 = 0°.

3 Results

3.1 Potential energy surface

The analysis of the fitted SAPT(DFT) PES reveals two minima. In Fig. 2 we show the global minimum geometry of the CS–H₂O SAPT(DFT) PES. It corresponds to the configuration in which the carbon atom of the CS molecule points toward one of the hydrogen atoms of the H₂O molecule, as is shown in Fig. 2, and corresponding to the minimum on curve C in Fig. 1. The c.o.m. separation is 8.11a₀, the angular coordinates are θ₁ = 122°, θ₂ = 297°, and ϕ₁ = ϕ₂ = 0°, and the well depth is 1219.23 cm⁻¹. A secondary minimum (curve marked as A in Fig. 1) is found for the geometry where the sulfur atom points toward the oxygen atom of the water molecule, with a c.o.m. separation of 7.06a₀ with θ₁ = θ₂ = ϕ₁ = ϕ₂ = 0° and a well depth of 660.44 cm⁻¹. The radial minimum for the arrangement B in Fig. 1 corresponds to the CS perpendicular to the H₂O plane, with the sulfur atom pointing toward H₂O. The c.o.m. separation for this minimum is R=7.36a₀, with angular coordinates of θ₁ = θ₂ = ϕ₁ = ϕ₂ = 90° and a well depth of 335.44 cm⁻¹. In Fig. 3, we show the anisotropy of the SAPT(DFT) PES for the cut showing potential as a function of θ₁ and θ₂ angular coordinates for a fixed R=8.11a₀, with other angles fixed at ϕ₁ = ϕ₂ = 0°.

It can be observed that the global minima and secondary minimum are separated by a barrier at about −328 cm⁻¹. The PES is found to be highly anisotropic with respect to both the θ₁ and θ₂ angles. The potential can change by about 1000 cm⁻¹ for a rotation of θ₁ smaller than 90°, and by roughly 1400 cm⁻¹ for a 180° rotation of θ₂. This shows that the interaction is very sensitive to the geometry of the complex.

Fig. 3 SAPT(DFT) PES of the CS–H2O complex for varying orientation angles, θ1 and θ2, with ϕ1 = ϕ2 = 0° at R = 8.11a0.

3.2 Rate coefficients

The behavior of the thermalized rate coefficients for CS–p-H₂O and CS–o-H₂O, obtained from Eq. (8) and later used in non-LTE excitation models, is shown in Fig. 4. The thermalized de-excitation rate coefficients generally decrease with increasing temperature, with a few exceptions such as the $j_{\mathrm{C}\mathrm{S}}=4\to j_{\text{CS\hspace{0.17em}}}^{\mathrm{\prime }}=0$ transition exhibited in Fig. 4, where the rate coefficient slightly increases from 5 to 20 K and then decreases smoothly up to 100 K. For both p-H₂O and o-H₂O, the favored transitions are the ones with the minimum kinetic energy release. In other words, as $\left|\mathrm{\Delta }{j}_{\mathrm{C}\mathrm{S}}\right|=|j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}-j_{\mathrm{C}\mathrm{S}}|$ increases, the magnitudes of the rate coefficients are decreasing. The CS–p-H₂O and CS–o-H₂O rate coefficients exhibit very similar behavior. In a systematic comparison (see Fig. B.1), CS–p-H₂O and CS–o-H₂O rate coefficients become increasingly similar as the temperature increases. At 10 K, deviations by more than a factor of two are observed for a few transitions, whereas at 100 K, the two sets of data agree within a factor of 1.5. In previous studies on the collisional excitation of CO, HF, and HCN induced by water molecules (Faure et al. 2020; Loreau et al. 2022; Żółtowski et al. 2025), the authors already observed that the excitation rate coefficients induced by p- and o-H₂O were relatively similar.

3.3 Excitation of CS in comets

The new sets of CS–H₂O thermalized rate coefficients were employed in a non-LTE radiative transfer model to investigate non-LTE effects on CS level populations in water-dominated comae, following the methodology described in Sect. 2.3. We also examined the impact of two approximations commonly used in the literature to compensate for the lack of collisional data. Cordiner et al. (2023) used the collisional rate coefficients of CS–p-H₂ and CS–o-H₂, computed by Denis-Alpizar et al. (2018), as substitutes for the CS–p-H₂O and CS–o-H₂O rate coefficients. The second approximation was employed by Biver et al. (1999), based on the work of Bockelée-Morvan (1987) and Crovisier (1987), who derived the rate coefficients (hereafter noted B99) from a single global cross section, σ_c = 5 × 10⁻¹⁴, cm², adapted from the H₂O-H₂O pressure broadening experiment of Murphy & Boggs (1969). The temperature-dependent rate coefficients were computed as $k_c(T)={\sigma }_c⟨v(T)⟩n_{j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}}(T)$, with ⟨v(T)⟩ being the mean velocity at the kinetic temperature, T, of the medium, and $n_{j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}}(T)$ being the population of the final CS level, assuming a Boltzmann distribution. In their approach, they do not differentiate between ortho and para water, and the rate coefficients inferred from this approximation are independent of the initial level.

In Fig. 5, we represent the populations of CS energy levels involved in the most frequently observed transitions – namely, j = 3 → 2, 5 → 4, and 7 → 6 (Roth et al. 2021a; Cordiner et al. 2023; Biver et al. 2023, 2024) – as a function of the nucleocentric distance, r. Therefore, the deviations observed in the populations are expected to be the ones most significantly affecting the modeling of CS emission spectra from the coma.

According to Fig. 5, the population of CS energy levels as a function of nucleocentric distance can be divided into three categories. Region (a) corresponds to the LTE domain, spanning distances from a few to a few hundred kilometers. This region represents the inner coma, where the density is high and molecular excitation is dominated by collisions. As a result, the CS level populations follow a Boltzmann distribution. Region (b) corresponds to the non-LTE domain, extending from a few hundred kilometers to approximately 10⁵ km from the nucleus. As the distance from the nucleus increases, the density decreases, and the population of CS energy levels deviates significantly from LTE conditions. Region (c) defines the fluorescence domain, where the excitation of CS is primarily governed by solar radiation. This corresponds to the outer part of the coma, where the density is so low that collisions can be neglected, as they no longer significantly affect the population of CS energy levels.

For single-dish observations, the telescope beam sizes are typically larger than a few arc-seconds, corresponding to nucleocentric distances greater than 10³ km for a comet at a geocentric distance of 1 a.u. (Bodewits et al. 2022). Therefore, the beam will necessarily include regions beyond the LTE domain (which, according to this study, extends to approximately 10² km), meaning that non-LTE modeling is required to accurately derive CS production rates in comets.

In Fig. 5, we also show the population of CS energy levels when CS–H₂ rate coefficients are used instead of CS–H₂O coefficients. We observe that the CS level populations differ significantly in the mid-coma region – where non-LTE effects are prominent – depending on which set of rate coefficients is used. Differences of over a factor of two can be observed between the populations, which can have a substantial impact on the modeling of CS production rates. This is not surprising, as a systematic comparison of state-to-state rate coefficients shows that CS–H₂O rate coefficients are significantly different than those of CS–H₂, by up to three orders of magnitude for the weakest transitions (see Fig. C.1). With CS–H₂O rate coefficients, the population of CS energy levels thermalizes at lower densities (i.e., at larger nucleocentric distances) due to more efficient energy transfer through collisions. This is due to the larger dipole moment of H₂O compared to H₂, resulting in stronger interactions with the target molecule, and so larger rate coefficients (see Fig. C.1).

In Fig. 5, we also show the population of CS energy levels using the set of B99 rate coefficients. The differences observed between the populations evaluated with the accurate CS–H₂O rate coefficients (computed in this work) and those obtained using the B99 datasets are very significant, and they can be by a factor of three. Thus, the use of B99 datasets is expected to lead to strong inaccuracies in the determination of the physical conditions within the coma.

We observe that the populations thermalize at slightly higher nucleocentric distances (i.e., lower densities), suggesting that the efficiency of collisions is overestimated with this approach. This is confirmed by the comparison between B99 and CS–H₂O rate coefficients, exhibited in Fig. D.1, where B99 data are generally greater than CS–H₂O data, and some transitions can be overestimated by over two orders of magnitudes at 100 K. Therefore, the analysis based on either approximation (using CS–H₂ or B99 datasets) should be revised, as using the accurate sets of data could lead to significant revision of CS production rates in comets.

In the work of Faure et al. (2018) on the excitation of CO by H₂O, only small differences were found between the CO populations derived using the accurate CO–H₂O dataset and those obtained with the global cross-section approximation of Biver et al. (1999) at 100 K. On the other hand, in the study on the excitation of HF by H₂O by Loreau et al. (2022), this approximation led to differences of up to a factor of 3 between the populations of HF obtained with the HF–H₂O rate coefficients they computed and one derived from the global cross sections employed by Bockelée-Morvan et al. (2014). Therefore, it remains uncertain on which collisional systems this approximation is expected to perform well, and it thus should be used with care.

Fig. 4 Representation of de-excitation thermalized rate coefficients (in cm3 s−1) of the CS–p-H2O (solid lines) and CS–o-H2O (dashed lines) from the first rotationally excited states (j1 = 1, 2, 3, 4) to the ground rotational state (j = 0).

Fig. 5 Population of CS energy levels (2 ≤ jCS ≤ 7) as a function of the nucleocentric distance (in km) using CS–H2O thermalized rate coefficients computed in this work (solid lines) and CS–H2 thermalized rate coefficients (dash-dotted lines) at 100 K. The range of nucleocentric delimits the following zones: (a) inner-coma, LTE; (b) mid-coma non-LTE; (c) outer coma fluorescence. Odd and even jCS have been separated for the sake of clarity.

4 Discussion

4.1 The expected accuracy of the CS–H₂O rate coefficients

To properly assess the accuracy of the CS–H₂O rate coefficients provided in this work, a comparison with experimental or with quantum CC calculations should be performed. As far as the authors are aware, no experimental measurements have been carried out on the CS–H₂O system, either of cross sections using crossed-beam experiments, or of rate coefficients, using double-resonances techniques for example. Regarding the CC calculations, the CS–H₂O system has a very deep potential well around 1200 cm⁻¹, meaning that many channels are coupled in the well, and thus that many closed channels are necessary to converge CC cross sections (this is also true with the coupled-states approximation). Therefore, it was not possible to converge transitions for the CS–H₂O system using the CC approach.

Another way to estimate the accuracy of the CS–H₂O rate coefficients is to compare the performance of SACM-derived rate coefficients for similar systems, for which comparisons with experimental or CC rate coefficients are available. The most suitable comparison system is CO–H₂O, as it has the same number of valence electrons, and its PES shares many similarities with that of CS–H₂O, particularly in the geometry of the global minimum – indicating similar chemical behavior.

For the calculation of the CO–H₂O rate coefficients, Loreau et al. (2018a) computed rate coefficients with the SACM method over the temperature range 5–100 K. The PES exhibits a global minimum around 650 cm⁻¹ – approximately half as deep as the minimum of the CS–H₂O PES, which lies at ~1200 cm⁻¹. As a result, the lifetime of the CO–H₂O complex is expected to be shorter than that of CS–H₂O. In addition, the number of available channels is significantly higher for CS–H₂O, as CS is heavier than CO. Therefore, the CS–H₂O intermediate complex is expected to exhibit a longer lifetime and a distribution of exit channels that more closely follows statistical behavior. Consequently, the SACM approach is expected to yield greater accuracy when applied to CS–H₂O collisions than to CO–H₂O. Loreau et al. (2018a) evaluated the accuracy of SACM cross sections by comparing them with full CC cross sections, limited to a few partial waves, J = 0–3, due to computational constraints⁴. They estimated that the SACM rate coefficients were accurate to within a factor of 1.5. Therefore, at least this level of accuracy – if not better – can be expected for the CS–H₂O rate coefficients presented in this work.

4.2 In the absence of cometary data: alternative strategies

In contexts in which collisional rate coefficients are unavailable – as in cometary atmospheres, where only five molecules (CO, HF, H₂O, HCN, and now CS) can be studied using non-LTE models with appropriate sets of rate coefficients – it is important to assess the accuracy of alternative strategies used to address this limitation. The following strategies can be employed to compensate for the lack of available rate coefficients in excitation models:

1. Deriving rate coefficients based on a global cross sections derived from pressure-broadening experiments (Bockelée-Morvan 1987; Crovisier 1987).

2. Using interstellar rate coefficients (with H₂ as a collider) to replace cometary ones (with CO, CO₂, or H₂O colliders, e.g., Cordiner et al. 2022, 2023)⁵.

3. Using rate coefficients involving other molecules based on the assumption that the collisional systems are similar enough (e.g., Cernicharo et al. 2024).

The first approximation is specific to the study of cometary comae and is widely used in contexts in which almost no collisional data exist. However, its impact on the derived physical parameters remains unknown. The second approximation assumes that excitation by light colliders (such as H₂) is presumably similar to that by heavier colliders (e.g., H₂O or CO). The final approximation is based on the assumption that if two molecules are sufficiently similar, their rate coefficients should also be similar. This approximation has been previously employed for interstellar applications (e.g., OCS–H₂ for HNCS–H₂ in Cernicharo et al. 2024), where the selection is based on data availability and similarities in the rotational structure between the two species. The impact of these approximations must be discussed in order to provide guidelines for cases in which rate coefficients are not available.

It was shown in the previous section that using either CS–H₂ or B99 rate coefficients instead of CS–H₂O rate coefficients is expected to substantially impact the population of CS energy levels, and consequently affect the physical conditions derived from non-LTE modeling based on these data. In this section the impact of using CO–H₂O rate coefficients instead of CS–H₂O is discussed. Indeed, both molecules are chemically similar, so the substitution seems reasonable.

The CO–H₂O rate coefficients, computed by Loreau et al. (2018a) and Faure et al. (2020) based on the PES of Kalugina et al. (2018), were taken from the EMAA database (Faure et al. 2025). Figure 6 presents a systematic comparison of the two datasets.

To assess the impact of the deviations on the non-LTE models in which these data will be employed and where dominant transitions will have more statistical weight, the weighted mean error factor (hereafter denoted WMEF) was calculated using Eq. (10), with $k_i^{\mathrm{R}\mathrm{E}\mathrm{F}}\equiv k_i^{\mathrm{C}\mathrm{S}-{\mathrm{H}}_2\mathrm{O}}$, and $k_i\equiv k_i^{\mathrm{C}\mathrm{O}-{\mathrm{H}}_2\mathrm{O}}$. $$\mathrm{W}\mathrm{M}\mathrm{E}\mathrm{F}=\frac{\sum _i{k}_i^{\mathrm{R}\mathrm{E}\mathrm{F}}{r}_i}{\sum _i{k}_i^{\mathrm{R}\mathrm{E}\mathrm{F}}},$$(10) where $r_i=max(k_i^{\mathrm{R}\mathrm{E}\mathrm{F}}/k_i,k_i/k_i^{\mathrm{R}\mathrm{E}\mathrm{F}})\ge 1$.

The CS–H₂O and CO–H₂O sets of rate coefficients are similar for both o- and p-H₂O, and at both temperatures. The differences are systematically less than a factor of two (except one transition at 50 K with o-H₂O). The CS–H₂O rate coefficients are generally larger than those of CO–H₂O, indicating more efficient collisions involving CS. This can be attributed to the larger dipole moment of CS compared to CO, which favors the energy transfer during collisions.

The impact of using CO–H₂O rate coefficients instead of CS–H₂O rate coefficients to model the population of CS energy levels was evaluated using the approach described in Sect. 2.3. Figure 7 represents the evolution of the population of the first CS energy levels with the nucleocentric distance using CS–H₂O thermalized rate coefficients computed in this work, and CO–H₂O thermalized rate coefficients computed by Loreau et al. (2018a) and Faure et al. (2020). As the CO–H₂O rate coefficients cover only the first 11 rotational levels of CO (up to j_CO ≤ 10), the CS–H₂O rate coefficient set was truncated to include the same number of levels in the non-LTE modeling. Therefore, only the effect of the rate coefficients themselves is discussed here.

We observe in this comparison that the population of CS energy levels along the whole range of nucleocentric distances is very similar when CS–H₂O and CO–H₂O data are employed. This is explained by the similarities between the two sets of rate coefficients, due to their chemical similarities, which are illustrated in Fig. 6. In the non-LTE domain (b) of Fig. 7, the population is observed to deviate from the fluorescence domain at shorter nucleocentric distances when CS–H₂O data are used, compared to CO–H₂O data. This is due to the larger CS–H₂O rate coefficients compared to the CO–H₂O ones. Similarly, the LTE domain (a) is reached at slightly larger nucleocentric distances for the same reason. Using CO as a template for CS in comets (and vice versa) seems to be a reasonable approximation, or at least a better approximation than using CS–H₂ rate coefficients or B99 rate coefficients. It is, however, difficult to predict in advance which set of collisional rate coefficients would serve as a good proxy. Therefore, any robust interpretation of cometary observations should rely on the use of the correct sets of collisional rate coefficients.

The fact that using CO–H₂O rate coefficients is a better approximation than using CS–H₂ data suggests that using the correct projectile is a better approach than considering the right target molecule. This findings suggest that the collisional partner (H₂O in our case) has more impact on the magnitude and behavior of the collisional data than the target molecule (here, CS).

To further assess the reliability of these approximations and the uncertainties they introduce into the physical parameters derived from observational analysis, they should be evaluated using real observations of various molecules. This would provide a clearer understanding of their impact on the determination of physical parameters, such as molecular production rates.

Fig. 6 Systematic comparison between the rate coefficients red(in cubic centimeters per second) of the CS–H2O system (x-axis) with the CO–H2O system (y-axis), with p-H2O in the left panel and o-H2O in the right panel) at 50 K (red crosses) and 100 K (black circles).

Fig. 7 Population of CS energy levels (2 ≤ jCS ≤ 7) as a function of the nucleocentric distance (in kilometers) using CS–H2O thermalized rate coefficients computed in this work (solid lines) and CO–H2O thermalized rate coefficients computed by Loreau et al. (2018a) and Faure et al. (2020, dashed lines) at 100 K. The range of nucleocentric delimits the following zones: (a) inner-coma, LTE; (b) mid-coma non-LTE; (c) outer coma fluorescence.

5 Conclusions

In this work, we have studied the collisional excitation of CS induced by both p-H₂O and o-H₂O in the 5–100 K temperature range, relevant for the study of cometary comae. The study began with the determination of the PES characterizing the electronic interaction between the two colliders. This PES was computed using the SAPT(DFT) and aVQZ level of theory. The PES was found to be highly anisotropic and presents a deep global minima at approximately 1220 cm⁻¹. The resulting 5D PES was implemented in the hibridon software to compute the adiabatic states of the collisional system. A probability was assigned to each adiabatic state following the SACM approach to compute the CS–p-H₂O and CS–o-H₂O state-to-state cross sections.

State-to-state rate coefficients were derived for transitions involving the first 19 rotational levels of CS (0 < j_CS < 18), the first five levels of p-H₂O (0₀₀ < j_kakc < 2₂₀), and the first five levels of o-H₂O (1₀₁ < j_kakc < 3₀₃), over the 5–100 K temperature range. The accuracy of the rate coefficients was discussed, and they are estimated to be accurate at better than a factor of 1.5 based on the validation of CO–H₂O data by Loreau et al. (2018a). Thermalized rate coefficients were calculated and used in non-LTE models to investigate the evolution of CS populations as a function of the nucleocentric distance. It was found that CS excitation is subject to strong non-LTE effects in the mid-coma, where excitation induced by solar radiation competes with collisional excitation induced by water molecules.

In addition, strategies for addressing the lack of rate coefficients in cometary atmospheres were examined. This study demonstrates that CS–H₂ rate coefficients are not reliable substitutes for CS–H₂O ones, nor is the global cross section suggested by Biver et al. (1999) a suitable alternative. Using CO–H₂O rate coefficients was found to be the best alternative. This suggests that selecting rate coefficients involving H₂O as the collisional partner and a chemically similar target molecule may yield more accurate modeling of the molecular excitation conditions than using data with the correct target molecule but with H₂ as the collisional partner. Further studies on other cometary systems are needed to confirm this finding, although such analyses are currently hindered by the limited availability of relevant data for cometary applications.

Data availability

The state-to-state rate coefficients produced in this work will be available on the BASECOL database [https://basecol.vamdc.eu, Dubernet et al. (2024)], and thermalized rate coefficients will be available on the LAMDA [https://home.strw.leidenuniv.nl/~moldata/, Schöier et al. (2005)] and EMAA [https://emaa.osug.fr, Faure et al. (2025)] databases.

Acknowledgements

AGP and FL thanks financial support from the European Research Council (Consolidator Grant COLLEXISM, Grant Agreement No. 811363), Rennes Metropole, and the PNP program (INSU CNRS). JK thanks for the financial support from the US National Science Foundation subaward No. 362627Sub1 of Grant No. AST-2009253. JK also is greatful for computational time provided on University of Maryland High-Performance Computing clusters. DK would like to thank the COST Action CA21101 “Confined Electronic Systems” (COSY) and gratefully acknowledges support from the National Science Centre (NCN) under Sonata Bis 9 grant No. 2019/34/E/ST4/00407. FL thanks the Institut Universitaire de France. The scattering calculations were performed with GENCI (project No. A0170413001), and the cluster of the Institut de Physique de Rennes (IPR).

References

Appendix A Computational details about the SACM method

In the SACM method, the rotational basis for the adiabatic states calculations should include all rotational levels of both monomers that are energetically accessible up to the maximum total energy considered for the collision. In our case, the goal was to compute rate coefficients between all j_CS ≤ 18 levels (at approximately 280 cm⁻¹), and all j_kakc ≤ 2₂₀ for p-H₂O, and up to j_kakc ≤ 3₀₃ for o-H₂O levels (both are lying at approximately 137 cm⁻¹). If we consider the kinetic energy distribution at 100 K to be of 400 cm⁻¹, then the highest total energy of the collision is the sum of these energies, so around 817 cm⁻¹.

According to the SACM method, any channels with an energy lower than this energy can be considered open, and have a non-zero probability in the adiabatic states counting. Therefore, any rotational level of both CS and H₂O with an energy lower than 817 cm⁻¹ should be included in the basis. Therefore, the ideal rotational basis would include all j_CS ≤ 36, and all j_kakc ≤ 10_kakc for both p-H₂O and o-H₂O. Note that the conversion between p- and o-H₂O being forbidden, the two sets of rate coefficients CS–p-H₂O and CS–o-H₂O are determined independently. However, employing this basis led to more than 10 000 channels for the partial wave J = 1, even with the EMAX option applied—surpassing the numerical limits typically acceptable for this type of scattering calculation. Consequently, truncation of the rotational basis was necessary.

The objective is therefore to reduce the size of the rotational basis without compromising the accuracy of the collisional rate coefficients obtained using the SACM approach. The truncated basis was finally set at j_CS ≤ 31, j_kakc ≤ 3_kakc for CS-p-H₂O, and j_CS ≤ 26, j_kakc ≤ 3_kakc for CS-o-H₂O. To ensure that this truncation was not affecting the quality of the data, the rate coefficients were computed with both rotational basis (full and truncated) for two partial wave as possible (J = 0, 1), and the two sets of rate coefficients are compared in Fig. A.1.

The WMEF was calculated based on Eq. (10), with $k_i^{\text{REF}}\equiv k_i^{\text{full}}$, the rate coefficients computed based on the adiabatic states computed with the full rotational basis, and $k_i\equiv k_i^{\text{trunc.}}$, the rate coefficients computed with the truncated rotational basis.

In this figure, we observe no significant differences between the two datasets, with a WMEF very close to unity in both cases. All dominant transitions are nearly perfectly reproduced by the truncated rotational basis, and the deviation remains within a factor of two at all times for the CS–o-H₂O rate coefficients. Only a few of the weakest transitions exceed a factor of two for CS–p-H₂O rate coefficients. Therefore, we conclude that this rotational basis can be reliably used to compute the CS–H₂O collisional data within the SACM approach.

Fig. A.1 Comparison of the (a) CS-p-H2O (b) CS-o-H2O rate coefficients (in cm3 s−1) computed with the full (on the x-axis) and truncated (on the y-axis) rotational basis for the adiabatic states calculations at 50 K (red and blue crosses) and 100 K (black circles).

Appendix B Comparison of CS–p-H₂O and CS–o-H₂O thermalized rate coefficients

Figure B.1 presents a systematic comparison between the CS–p-H₂O and CS–o-H₂O thermalized rate coefficients at 10 K and 100 K. The WMEF was calculated based on Eq. (10), with $k_i^{\mathrm{R}\mathrm{E}\mathrm{F}}\equiv k_i^{\mathrm{C}\mathrm{S}-p-{\mathrm{H}}_2\mathrm{O}}$, the CS–p-H₂O rate coefficients, and $k_i\equiv k_i^{\mathrm{C}\mathrm{S}-o-{\mathrm{H}}_2\mathrm{O}}$, the CS–o-H₂O rate coefficients.

It can be observed that both sets are relatively similar, with agreement always better than a factor of two at 100 K, and with nearly all transitions agreeing within a factor of two at 10 K. The WMEF values—1.36 at 10 K and 1.09 at 100 K—indicate that the dominant transitions are similar in both cases. The agreement between the o- and p-H₂O rate coefficients improves as the temperature increases, which can be attributed to a reduced impact of the energy threshold as the temperature increases. This effect has also been observed and discussed by Żółtowski et al. (2025).

With regard to the similarities between the two datasets, it is unlikely that the ortho-to-para ratio chosen in the non-LTE modeling will significantly impact the modeling of CS in water-dominated comets.

Fig. B.1 Systematic comparison between thermalized de-excitation rate coefficients of CS–p-H2O and CS–o-H2O collisional system at 10 K (green crosses) and 100 K (black circles).

Appendix C Comparison of CS-H₂O and CS-H₂ rate coefficients

For the modeling of CS in water-dominated comets, the collisional rate coefficients of CS with H₂O—previously unavailable—have sometimes been substituted with the CS–H₂ rate coefficients computed by Denis-Alpizar et al. (2018) using the CC approach. To assess the accuracy of this approximation, the CS–H₂O rate coefficients computed in this work are compared to the CS–H₂ rate coefficients taken from the EMAA database (Faure et al. 2025) in Fig. C.1.

The deviations observed between the two sets of rate coefficients—those involving H₂O and H₂ as projectiles—are very significant. The dominant transitions tend to be stronger for CS–H₂ compared to CS–H₂O for both ortho- and para-H₂O, whereas the weaker transitions are generally much stronger for CS-H₂O thermalized rate coefficients. This is consistent with the behavior of rate coefficients involving light projectiles such as H₂, which tend to decrease rapidly with increasing Δj_CS, while remaining relatively stable for heavier colliders like o/p-H₂O.

The PES of the two collisional systems differ significantly in their well depths—the well for the CS–H₂O system is much deeper than that for CS–H₂ (1269 cm⁻¹ vs. 173 cm⁻¹). Since the intermediate complex formed by CS and H₂O has a longer lifetime than the one formed with H₂, the efficiency of the collision becomes less sensitive to the amount of energy transferred, and therefore less dependent on Δj_CS.

At 20 K, the global deviation, as indicated by the WMEF, is about 13.7 for p-H₂O vs. p-H₂ and 8.07 for o-H₂O vs. o-H₂. These deviations decrease to approximately a factor of 4 for both o/p-H₂O vs. o/p-H₂ at 100 K. This improvement can be attributed to the fact that, at higher temperatures, the rate coefficients involving H₂ as a projectile tend to vary less with quantum numbers, thus aligning more closely with the behavior of the CS–H₂O rate coefficients and reducing the overall deviation.

However, it is clear that the CS–H₂ rate coefficients are not a reliable substitute for the CS–H₂O rate coefficients. Using them in place of the appropriate data is likely to introduce significant inaccuracies in the derived physical conditions when modeling CS emission spectra in cometary comae.

Fig. C.1 Systematic comparison between the rate coefficients (in cm3 s−1) of CS in collision with H2O (x axes) and H2 (y-axis) (p–H2/p–H2O in the left panel and o–H2/o–H2O in the right panel) at 20 K (blue crosses) and 100 K (black circles).

Appendix D Comparison of CS-H₂O and B99 rate coefficients

The global cross-section approach has been proposed by Bockelée-Morvan (1987) and Crovisier (1987), who derived the rate coefficients from a single global cross section adapted from the H₂O-H₂O pressure broadening experiment of Murphy & Boggs (1969). For CS–H₂O rate coefficients, Biver et al. (1999) using a global cross section of σ_c = 5 × 10⁻¹⁴, cm², from which they obtained B99 rate coefficients as $k_c(T)={\sigma }_c⟨v(T)⟩n_{j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}}(T)$, with ⟨v(T)⟩ being the mean velocity at the kinetic temperature T of the medium, and $n_{j_{\mathrm{C}\mathrm{S}}^{\mathrm{\prime }}}(T)$ being the population of the final CS level, assuming a Boltzmann distribution. In their approach, they do not differentiate between ortho and para water, and the rate coefficients inferred from this approximation are independent of the initial level.

The B99 rate coefficients are compared to the thermalized CS–H₂O rate coefficients obtained in this work in Fig. D.1.

The deviation between the two sets of rate coefficients are particularly large, especially at 20 K where the WMEF indicate over two order of magnitude between the sets of data. This is due to the thermal weight attributed, which is going to be very low for certain CS level at 20 K, but remain relatively equally spread among CS energy levels at 100 K, thus leading to a better agreement with CS–H₂O data.

It is important to note that B99 rate coefficients are systematically higher than the CS–H₂O rate coefficients, suggesting an overestimation of the global cross section proposed by Biver et al. (1999). Such differences are still very significant, with a WMEF equal to 4.21, suggesting that this approximation is not expected to perform well to accurately constrain physical parameters from CS molecular spectra in comets.

Fig. D.1 Systematic comparison between rate coefficients (in cm3 s−1) of CS in collision with thermalized H2O (x axes) and B99 rate coefficients (y-axis) at 20 K (magenta crosses) and 100 K (black circles).

All Figures

	Fig. 1 Comparison of the CS–H2O SAPT(DFT) and aVQZ PES fit (solid black line) with the CCSD(T)-F12 and aVTZ+BF calculations performed for selected radial cuts. The θ1 and ϕ1 are spherical angles of the R vector, while θ2 and ϕ2 are spherical angles describing the orientation of the CS molecule.
In the text

	Fig. 2 Global minimum equilibrium geometry of the CS–H2O SAPT(DFT) and aVQZ PES. The Re between the c.o.m. of H2O and c.o.m. of CS is 8.11a0, with θ1 = 122°, θ2 = 297°, and ϕ1 = ϕ2 = 0°.
In the text

	Fig. 3 SAPT(DFT) PES of the CS–H2O complex for varying orientation angles, θ1 and θ2, with ϕ1 = ϕ2 = 0° at R = 8.11a0.
In the text

	Fig. 4 Representation of de-excitation thermalized rate coefficients (in cm3 s−1) of the CS–p-H2O (solid lines) and CS–o-H2O (dashed lines) from the first rotationally excited states (j1 = 1, 2, 3, 4) to the ground rotational state (j = 0).
In the text

Fig. 5 Population of CS energy levels (2 ≤ jCS ≤ 7) as a function of the nucleocentric distance (in km) using CS–H2O thermalized rate coefficients computed in this work (solid lines) and CS–H2 thermalized rate coefficients (dash-dotted lines) at 100 K. The range of nucleocentric delimits the following zones: (a) inner-coma, LTE; (b) mid-coma non-LTE; (c) outer coma fluorescence. Odd and even jCS have been separated for the sake of clarity.

In the text

	Fig. 6 Systematic comparison between the rate coefficients red(in cubic centimeters per second) of the CS–H2O system (x-axis) with the CO–H2O system (y-axis), with p-H2O in the left panel and o-H2O in the right panel) at 50 K (red crosses) and 100 K (black circles).
In the text

Fig. 7 Population of CS energy levels (2 ≤ jCS ≤ 7) as a function of the nucleocentric distance (in kilometers) using CS–H2O thermalized rate coefficients computed in this work (solid lines) and CO–H2O thermalized rate coefficients computed by Loreau et al. (2018a) and Faure et al. (2020, dashed lines) at 100 K. The range of nucleocentric delimits the following zones: (a) inner-coma, LTE; (b) mid-coma non-LTE; (c) outer coma fluorescence.

In the text

	Fig. A.1 Comparison of the (a) CS-p-H2O (b) CS-o-H2O rate coefficients (in cm3 s−1) computed with the full (on the x-axis) and truncated (on the y-axis) rotational basis for the adiabatic states calculations at 50 K (red and blue crosses) and 100 K (black circles).
In the text

	Fig. B.1 Systematic comparison between thermalized de-excitation rate coefficients of CS–p-H2O and CS–o-H2O collisional system at 10 K (green crosses) and 100 K (black circles).
In the text

	Fig. C.1 Systematic comparison between the rate coefficients (in cm3 s−1) of CS in collision with H2O (x axes) and H2 (y-axis) (p–H2/p–H2O in the left panel and o–H2/o–H2O in the right panel) at 20 K (blue crosses) and 100 K (black circles).
In the text

	Fig. D.1 Systematic comparison between rate coefficients (in cm3 s−1) of CS in collision with thermalized H2O (x axes) and B99 rate coefficients (y-axis) at 20 K (magenta crosses) and 100 K (black circles).
In the text