© The Authors 2025

1 Introduction

The first detection of diatomic carbon by Souza & Lutz (1977) towards Cygnus OB2-12 was obtained in the near infrared (1-0) Phillips band. That study was followed by several other detections in the same electronic system. Lambert et al. (1995) made definitive interstellar detections of the D - X and F - X electronic transitions from ultraviolet observations with the Hubble Space Telescope. Sonnentrucker et al. (2007) and Hupe et al. (2012) expanded the sample of sight lines with UV measurements. These observations allow us to infer the column density of the different rotational levels up to relatively high values of J and are thus promising diagnostic tools of the physical conditions of the associated environments if the different excitation mechanisms are well accounted for. Since then, higher resolution observations have been obtained in the visible and near infrared with the Akari telescope (Hamano et al. 2019), the very large telescope ultraviolet and visual echelle spectrograph (VLT UVES) (Fan et al. 2024), and the new habitable zone planet finder (HPF) (present work). The Hobby-Eberly telescope (HET) HPF spectra significantly extend observations of the v = 0 X state rotational ladder beyond previous upper limits, which opens the possibility of novel tests.

van Dishoeck & Black (1982) presented the first detailed analysis of C₂ excitation by computing the quadrupole radiative transition probabilities within the ground state up to J = 16, and the intercombination transitions between levels of the lowest triplet and singlet states, as well as estimating the rotational collisional de-excitation rates due to collisions with H and H₂. Le Bourlot & Roueff (1986) improved the treatment of the intercombination transitions and introduced them into a model applied to the ζ Oph line of sight.

More recently, Casu & Cecchi-Pestellini (2012) included new spectroscopic results on the C₂ Phillips band and collisional excitation rates due to He, and proposed a two-phase model that enabled them to explain the observations in a few diffuse and translucent lines of sight. Finally, Neufeld et al. (2024) introduce the new collisional rate coefficients for the excitation of C₂ due to para- and ortho-molecular hydrogen (Najar & Kalugina 2020) by using the spectroscopic information reported in Casu & Cecchi-Pestellini (2012). Their study leads to a significant reduction in the density estimates as a consequence of the larger collisional rate coefficients computed by Najar & Kalugina (2020).

In this paper, we revisit the excitation mechanisms pertaining to interstellar¹ C₂ and present a 0D excitation model that is further compared to new observations. In Sect. 2 we describe the detailed balance equations and the different processes that contribute to the population of C₂ rotational levels, including collisional, radiative, and chemical formation contributions. Sect. 3.1 details the importance of the various contributions for a given set of temperatures and densities. In Sect. 4 we report the observations of C₂ in Cyg OB2-12, HD 29647 performed with the HPF and in HD 24534 (Sonnentrucker et al. 2007) and in Sect. 5 we present our derivation of physical conditions through an optimization schema of our excitation model. Our conclusions are summarized in Sect. 6.

Fig. 1 Top: C2 electronic systems included in this computation. “I.C.” stands for “Intercombination” (inspired by Yurchenko et al. 2018). Bottom: Ro-vibrational levels included in the full detailed balance computation. Numerical values are given in Table A.1.

2 Excitation mechanisms of C₂

2.1 Detailed balance equations

It is customary to study the detailed balance of a species by only considering collisions and radiative transitions within the set of levels considered. However, this leads to serious inaccuracies in the resulting populations if other physical processes occur with characteristic times comparable to collision times or radiative lifetimes. In the case of H₃⁺, Le Bourlot et al. (2024) have shown that chemical destruction by electrons may be fast enough to prevent thermalization of the lowest ortho level of the molecule, leading to an excitation temperature between the two lowest levels that is typically half that of the gas kinetic temperature.

These state-specific chemical rates are easily included in the balance equation and have the advantage, for the theoretician, of providing a solution without using mass conservation. Hence we have the further benefit of a direct check of the quality of the solution. The general form of the balance equation for a level i of any species is thus $$x_i(\sum _j{R}_{ij}+D_i)=\sum _j{x}_j{R}_{ji}+F_i$$(1)

Here, x_i is the relative abundance of level i, R_ij is a rate of transition (including inelastic collisions and radiative decay or pumping) from level i to level j (respectively from j to i) in s⁻¹, D_i is the destruction rate of C₂ on level i by “chemical” processes in s⁻¹, and F_i is the formation rate of level i in s⁻¹. It is important to note that excitation can never be decoupled from other physical processes. State-specific formation and destruction must always be accounted for. These processes are included here for the first time in C₂.

The following subsections detail the different contributions to R_ij. The electronic $X{}^1{\mathrm{\Sigma }}_g^{+}$ ground state of C₂ is very close to the a ³Πu excited state so that to compute accurate ground-state populations of rotational levels up to J = 32, as found towards some lines of sight, one must also include rotational levels belonging to a ³Πu that fall between rotational levels of X. In addition, vibrationally excited levels with a low rotational quantum number belonging to X may decay towards a². Finally, we computed in a consistent way the detailed balance of the first 120 levels that encompass the v = 0, J = 34 level of X, in order to avoid border effects.

The top panel of Fig. 1 shows the electronic states included in this computation through radiative pumping followed by cascades (see Sect. 2.4); the bottom panels shows the position of all 120 levels included in the full computation. This includes 18 rotational levels within X, v = 0, 8 rotational levels within X, v = 1, and 94 levels within a. One can see that the density of levels is much higher in the a state, which leads to a dilution of abundances.

2.2 Collisions

The C₂ molecule forms mainly in the molecular part of clouds. Hence, the main collision partner by far is H₂. Najar & Kalugina (2020) provide separate rate coefficients for para- and ortho-H₂ collisions within X for J  20. This is not enough to compute accurate populations up to J = 32. Thus, we have extrapolated these rates for Δ J = 2, 4, 6, and 8 with a simple exponential fit. Although rather crude, this approximation is better than setting the rates to 0. It is shown on Fig. 2 and numerical values are provided in Tables A.4 and A.5.

Collisional excitation rate coefficients due to He are found in Najar et al. (2008) (Table 1) for T = 5 to 100 K upto J = 20 and in Ritika & Dhilip Kumar (2022) for very low temperatures and a few J. They agree roughly at 5 K. The same procedure as for collisional rates with H₂ was used to extend the available rates to higher J for Δ J = 2 and 4, as found in Table A.6.

To the best of our knowledge, no collision rates exist for collisions with H, so all models below are computed for a fully molecular gas. If collision rates with H are similar to those with para-H₂ (as is often the case), then a higher density n_H would be required to induce the same impact of collisions on excitation. For collision-induced transitions, $R_{ij}=k_{ij}^{\mathrm{X}}{n}_{\mathrm{X}}$, where $k_{ij}^{\mathrm{X}}$ is the rate coefficient for collisions with species X of abundance n_X.

Fig. 2 Extrapolation of collisional rate coefficients from Najar et al. (2008) at 30K. Extrapolated rate coefficients are on the right side of the black vertical line. Top: Para-H2, Bottom: Ortho-H2.

2.3 Electric quadrupolar radiative transitions within X ground electronic state

As C₂ is a homonuclear molecule, there are no electric dipolar transitions within the X electronic state. Quadrupolar transition probabilities have been computed by van Dishoeck & Black (1982) but only up to J = 20 (their Table 2). Here, we propose to use the energy dependence of these expressions to justify an extrapolation scheme. Electric quadrupole emission coefficients for J′ → J″ are given by $$A_{J^{\prime },J^{″}}=\frac{32{\pi }^6}{5h}{\overline{\nu }}^5\frac{S_Q(J^{\prime },J^{″})}{2J^{\prime }+1}$$(2)

where ν̄ is the transition wavenumber in cm ⁻¹ and the strength S_Q (J′, J″) is in atomic units: $$S_Q(J^{\prime },J^{″})={|\int {\psi }_{J^{\prime }}^{\ast }Q{\psi }_{J^{″}}dR|}^2\frac{J^{\prime }(J^{\prime }-1)}{2J^{\prime }-1}$$(3)

$\overline{\nu }(J^{\prime }\to J^{\prime }-2)=4B{J}^{\prime }-2B$, with B being the rotational constant of the molecule, so that ${\overline{\nu }}^5\propto J^5$. From Eq. (3) we see that, for high rotational J values, $\frac{S_Q(J^{\prime },J^{″})}{2J^{\prime }+1}$ does not depend on J. So A_J′,J″ varies as J⁵ (neglecting the possible variation of the electric quadrupolar transition matrix element). Plotting log₁₀ (A_J′,J″) as a function of log₁₀(J), we see that the points corresponding to the highest J transitions are aligned on a straight line (Fig. 3). The transition probability values extrapolated from a linear fit and the transition wavenumbers obtained from the ExoMol database³ are reported in Table A.2.

For v > 0 transitions, van Dishoeck & Black (1982) only give the sum A (v′, J″) over all lower levels, and state that these A (v, J) are “virtually independent of J”. Using the expressions of S (J′, J″) for the S, Q, and O branches, we can distribute these sums to recover individual A_{v′=1J′,v″=0J″}, as displayed in Table A.3. For electric quadrupolar transitions, R_ij = A_j, the Einstein emission coefficient.

Fig. 3 Electric quadrupole radiative rotational transition probabilities from van Dishoeck & Black (1982). The positions of J = 22 and 24 are shown by straight lines. The points above are extrapolations.

2.4 Electronic radiative processes

In their Fig. 2, van Dishoeck & Black (1982) describe the different radiative electronic transitions that link the rotational levels of the $X{}^1{\mathrm{\Sigma }}_g^{+}$ ground electronic state. However, in the absence of complete radiative information, these authors introduced some approximations to account for these mechanisms. Radiative electronic transitions data within the eight lowest electronic states of C₂ are now accessible and reported in the ExoMol database, extracted from the recent study of McKemmish et al. (2020). The database includes transitions within singlet states $X{}^1{\mathrm{\Sigma }}_g^{+},A{}^1{\mathrm{\Pi }}_u,B{}^1{\mathrm{\Delta }}_g$, and $B^{\prime }{}^1{\mathrm{\Sigma }}_g^{+}$ and within triplet states a³Πu, $b{}^3{\mathrm{\Sigma }}_g^{+}$, $c{}^3{\mathrm{\Sigma }}_u^{-}$, and d³Π_g. The ExoMol database also incorporates intercombination transitions between singlet and triplet states. We have supplemented it with transitions to $D{}^1{\mathrm{\Sigma }}_u^{+}$ (Mulliken system) with data from Schmidt & Bacskay (2007) (oscillator strengths) and Krechkivska et al. (2018) (energy level positions and derived transition wavenumbers). This leads to a set of over 40 000 levels coupled by more than 600 000 lines. Within this total of N = 44 189 levels, we only computed a detailed balance for the lowest M ones (M = 120 here to include all levels up to X, v = 0, J = 34, just above the highest level observed; see Sect. 5.2).

These data provide the spontaneous radiative emission coefficients A_ij and corresponding transition wavenumbers v̄, from which absorption B_ji and induced emission B_ij coefficients can be derived. In the presence of an external radiation field, such as the ambient interstellar radiation field (ISRF), absorption toward electronically excited states is possible, followed by radiative cascades to a low lying level. This mechanism results in the possible coupling of two low lying levels even if no direct transition is possible between them. For two levels i and j below M, the coupling coefficient by this indirect cascade, R_ij, can be written in the optically thin approximation as $$R_{ij}=\sum _{n>M}{B}_{in}{\overline{J}}_{in}{C}_{nj}$$(4)

where B_in is the Einstein absorption coefficient to a level n above M, J̄_in is the mean radiation field intensity, and C_nj is the cascade coefficient from n to j. For a fixed value of M, the coefficients C_nj may be computed in advance and stored. With M = 120, we end up with N = 12 878 levels, which are accessible through 71 708 radiative transitions and lead to 1 463 887 cascade coefficients. The algorithm to compute the cascade coefficients is given in Appendix B.

Fig. 4 Adopted ISRF for G = 1 in the range of possible transitions of C2. Colored points displayed on the horizontal line at the bottom of the figure indicate the positions of included transitions coded with the value of Blu Jul from red (strongest) to light green (weakest).

2.5 Adopted ISRF

The adopted ISRF is the one used in the Meudon PDR code⁴, based on Mathis et al. (1983). It is defined as a reference and corresponds to G=1. That dimensionless factor will be used further for possible future scalings of the ISRF. However, the scaling does not apply in the millimeter domain where the Cosmic Microwave Background (CMB) contribution applies. It is displayed in Fig. 4. It extends from the ultraviolet (where transitions from the Mulliken system take place) to the submillimeter range for the lowest rotational transitions. Neufeld et al. (2024) only consider the 0.77-1.21 μm restricted spectral window, covering the Phillips $A{}^1{\mathrm{\Pi }}_u-X{}^1{\mathrm{\Sigma }}_g$ band system of C₂, as shown on top of the present adopted curve in Fig. 4, where an almost perfect agreement is found with our adopted ISRF in the considered spectral range.

However, additional noticeable absorptions take place outside that window. Horizontal points on Fig. 4 show the positions of these absorptions and are color-coded proportionally to B_lu J_ul. In particular, weak but non-negligible absorptions also occurs in the mid-infrared region at the positions of some polycyclic aromatic hydrocarbons (PAH) features.

We computed the resulting absorption probability of the X and a levels included in our excitation model for G = 1 and found an almost J independent value of 4.3 10⁻⁹ s⁻¹ for the rotational X levels and obtained a value of ~6.7 10⁻⁹ s⁻¹ for the a levels. We note that Neufeld et al. (2024) mention this same rotational independent value for their estimated pumping rate of the X rotational levels that results from their update of the radiation field used previously by van Dishoeck & Black (1982).

2.6 Chemical excitation

In the cold ISM, C₂ is mostly formed by exothermic reactions. Detailed examination of the various formation paths shows that no single reaction dominates, but that a few contribute at the 10 to 20% level. Among these, some reactions have a large exothermicity. The main ones are: $$\begin{array}{r}{{\mathrm{C}}_2\mathrm{H}}^{+}+e^{-}\to {\mathrm{C}}_2+\mathrm{H}(\mathrm{\Delta }H=6.671\mathrm{e}\mathrm{V}),\\ {{\mathrm{C}}_3}^{+}+e^{-}\to {\mathrm{C}}_2+\mathrm{C}(\mathrm{\Delta }H=4.804\mathrm{e}\mathrm{V}),\\ {\mathrm{C}}^{+}+\mathrm{C}\mathrm{H}\to {\mathrm{C}}_2+{\mathrm{H}}^{+}(\mathrm{\Delta }H=0.321\mathrm{e}\mathrm{V}).\end{array}$$

Electronic recombination rates are found in Florescu-Mitchell & Mitchell (2006). Branching ratios for light carbonaceous molecules are found in Chabot et al. (2013). Level number 120 of C₂, corresponding to F1, a ³Πu, v = 1, J = 2, is about 3000 K above the ground level. So, if a significant fraction of the exother-micity remains in the molecule, it is enough to populate levels proportionally to their statistical weight. This is a zeroth order approximation, which allows for easy testing of the importance of chemical formation excitation, compatible with the current state of knowledge⁵.

In the following, we take $$F_i=k_f\frac{g_i}{\sum _M{g}_j}$$(5)

where k_f (in s⁻¹) is a single effective parameter that accounts for all chemical reactions leading to internal excitation and g_i is the statistical weight of level i (f is for “formation”). All M levels are included, both from X and a. We define k_f by $$k_f=\frac{1}{[{\mathrm{C}}_2]}{\frac{d[{\mathrm{C}}_2]}{dt}|}_f$$(6)

Optimal values for k_f derived below have been checked for consistency with the results from the Meudon PDR code, where the full chemistry is solved, including over 8000 reactions between over 230 species. Typical values go from a few 10⁻¹³ to 10⁻¹¹ s⁻¹, depending on cloud structures and conditions.

No simple reasoning leads to a particular choice of statespecific destruction rates. So we take simply D_i = k_f (in s⁻¹), which ensures that $$\sum _i{x}_i{D}_i=\sum _i{F}_i.$$(7)

3 Results from 0D model

We now discuss the numerical results⁶ obtained after including the various excitation processes in the balance equations (Eq. (1)). The input parameters are the density, n_H, the temperature, T, the ambient radiation field scaling, G, and the chemical formation rate, k_f. The output values are the relative populations of C₂.

The system of M equations is linear and well behaved, and no conservation equation is needed if k_f ≠ 0. It is easily solved using, e.g., the LAPACK library. However, using 120 levels leads to a system with a rather high condition number and requires the use of the DGESVX extended routine to allow for matrix equilibration. We suppose that hydrogen is fully molecular, so that n (H₂) = n_H/2, and that the ortho to para ratio of H₂ is at thermal equilibrium. This assumption is fine for diffuse cloud conditions where the gas temperature is derived from the J=1 to J=0 ratio of H₂ (Shull et al. 2021; Le Bourlot et al. 2024). We also include helium, which represents 20% of the molecular H₂ density. Here we emphasize the very large dependency of the collisional excitation rate coefficients of C₂ on the ortho to para ratio of H₂ , as already pointed out by Neufeld et al. (2024). Such a 0D model is restricted to study the relative populations of a defined number of C₂ levels without considering the detailed chemical or thermal balance and independent of geometry. It may be considered as an extension and update to the Ben McCall C₂ Calculator that includes only seven levels of C₂⁷.

Fig. 5 C2 excitation diagram for T = 30K and nH = 50 cm−3 without radiative cascades. Red: Only levels from X included (labeled “X only”); Blue: X and a electronic states, including intercombination transitions (labeled “X + a”).

3.1 A fiducial example

The impact of these various excitation mechanisms is best seen on a representative example. We consider gas at a temperature of T = 30K and n_H = 50 cm⁻³ and introduce the different excitation mechanisms step by step. In the following, we assume a radiation field scaled by the value G = 1. We first examine the dependence of the weighted fractional abundances of C₂ on the inclusion of state a. Only collisions and direct radiative absorption or emission are accounted for (Fig. 5). Using only the rotational levels from X (red points), we recover a pure Boltzmann distribution at the gas temperature. This is because quadrupolar transitions are weak and even a density of 50 cm⁻³ is above the critical density. Adding the a levels, which are in the same energy range as the rotational levels belonging to X (the previously mentioned 120 levels), leads to the weighted populations in blue. Radiative pumping from low-lying rotational levels of X gives a perceptible population of state a showing up at an energy of ~800 K. The fractional populations within X remain the same.

Fig. 6 is restricted to the X rotational levels, although levels from a are included in the computation. This shows the effect of adding further excitation mechanisms. Blue points (label X) are the same as in Fig. 5. Green points show the effect of radiative pumping from Sect. 2.4 (label “X + Cascade”). This approximation represents a significant update from the previous treatments (van Dishoeck & Black 1982; Casu & Cecchi-Pestellini 2012; Neufeld et al. 2024), where intercombination transitions are included approximately over a restricted spectral range. The corresponding weighted fractional abundances of C₂ decrease smoothly as a function of the rotational energy. We have also tested this case with and without the inclusion of the Mulliken electronic band system in the cascade coefficients. The differences are insignificant up to J = 30 and increase to only 10% for J= 34.

Light blue points refer to a computation without considering radiative cascades but with chemical pumping (Sect. 2.6) with k_f = 8 × 10⁻¹² s⁻¹ (label “X + Chem”). Chemical pumping maintains an almost constant abundance for J > 6. Red points include all contributions together (label “X + Cascade + Chem”). Radiative cascades dominate for intermediate J (from J = 6 to 12); then chemical excitation maintains a slowly decreasing population.

Fig. 6 Effect of nonthermal excitation mechanisms on the weighted fractional populations (see text for the explanation of the labels).

Fig. 7 Excitation diagram for a ratio nH/G = 50. All results collapse on a single curve up to J = 16. Note that the density range is much larger than in the cases of Appendix C.

3.2 Parameter dependencies

Such a model is hopefully useful to extract the relevant physical conditions from the excitation properties. In the following, we discuss the dependence of these on density, temperature, and radiation field. We consider chemical pumping effects that impact only the most rotationally excited levels separately.

Fig. 7 shows an excitation diagram computed for T = 30K, k_f = 0 and a range of densities n_H. However, we have set a fixed value to the ratio α = n_H/G. We see that all curves collapse on a single excitation diagram up to level J = 16. This may be understood by considering the various terms of Eq. (1). We observe that quadrupolar spontaneous transitions A_ij are much smaller than other terms for J < 16, and can be neglected. If we additionally neglect chemical terms D_i and F_i, then Eq. (1) reduces to $$n_i\sum _j(k_{ij}^X{n}_X+\sum _{n>M}{B}_{in}{\overline{J}}_{in}{C}_{nj})=\sum _j{n}_j(k_{ji}^X{n}_X+\sum _{n>M}{B}_{jn}{\overline{J}}_{jn}{C}_{ni})$$(8)

The density of collision partners n_X is proportional to n_H, while the mean radiation field J̄ is proportional to G. Thus, substitution of n_H by α G allows for a cancellation by G. The resulting populations only depend on α. Overcoming this degeneracy requires additional information. The presence of nearby stars, in particular, may impact the strength and the shape of the radiation field. The resulting population departs from this scaling if one of the neglected terms starts to play a role. This is the case for high J levels, for which the emission coefficients A_ij are large enough to contribute, as seen in Fig. 7 (see Table A.2). Including chemical pumping (k_f ≠ 0) also lifts the degeneracy, since k_f does not scale with n_H (see Eq. (6)). As seen below, its effect on high J levels is much more efficient than radiative pumping and cascades alone. We have checked that this result remains true for all relevant temperatures (between 10 and 100 K).

In Appendix C, we study the sensitivity of the excitation diagram to systematic individual variations of the physical parameters, temperature, density, and strength of the radiation field. It is clearly shown that none of those enables us to obtain enhanced populations for J values larger than ~18, as found in the observations. Then, chemical excitation appears as the only mechanism that allows such a nearly constant ratio N_J/g_J, as observed for high Js.

4 Observations

We consider three different lines of sight in order to appraise our excitation model against observations. Observations toward HD 29647 and Cyg OB2 12 were acquired with the HPF on the 10 m HET at the McDonald Observatory (Mahadevan et al. 2014; Mahadevan et al. 2018) as part of a larger sample to examine C₂ and CN excitation and isotopologue ratios at near-infrared wavelengths. Here we provide a brief summary of the data acquisition and processing, leaving the detailed description to a future paper (Federman et al., in prep.). Besides these targets, spectra of telluric standards were obtained to remove lines arising from Earth’s atmosphere. The C₂ lines from the (2-0) and (1-0) bands of the Phillips system were used in the present analyses. The detector was a 2048 × 2048 pixel Hawaii-2 HgCdTe device with a resolving power of about 53 000. The software package Goldilocks⁸ was used for data reduction. The extracted onedimensional spectra were converted to air wavelengths through the use of Eq. (16) in Morton (2003). After dividing the target spectrum by that of the telluric standard, the target spectrum was normalized and placed on the local standard of rest. The final spectra had signal-to-noise ratios per resolution element of about 2500.

The C₂ spectra from HPF were fit using the profile synthesis code ismod (Sheffer et al. 2008), from which the velocity (v_LSR), column density for each rotational level (N(J)), and Doppler parameter (b-value) were obtained. We used the experimental measurements on C₂ wavelengths from Chauville et al. (1977) and the results of theoretical results for oscillator strength by Schmidt & Bacskay (2007). For the gas toward HD 29647 and Cyg OB2 12, as starting points we adopted the previously determined component structures by Federman et al. (2021) and Hamano et al. (2019), respectively. Our results for HD 29647 are generally consistent with earlier measurements (Federman et al. 2021), and the agreement between our results for Cyg OB 12 and those of Hamano et al. (2019) is very good (Federman et al., in prep.). Since the data for modeling C₂ excitation came from the ExoMol website based on recent calculations by McKemmish et al. (2020), updating the values from Yurchenko et al. (2018), we compared these results with the relevant values used in fitting the spectra to check for consistency. The wavelengths are consistent, considering the spectral resolution of HPF, and the line oscillator strengths agree at the 1 to 2% level.

However, we note another point regarding the transition wavenumbers and associated air wavelengths of Q lines displayed in the ExoMol database, where the discrepancies between experimental (Chauville et al. 1977) and reported values in Exo-Mol increase with the rotational quantum number. This is likely due to the presence of Λ-doubling, which was not taken into account in the ExoMol calculations (Tennyson, Yurchenko, priv. comm.).

The data for C₂ toward HD 24534 are taken from Table 1 in Neufeld et al. (2024) and correspond to the average value derived from observations in the ultraviolet F-X and D-X systems with the Space Telescope Imaging Spectrograph (STIS) and visible Phillips system with the Astrophysical Research Consortium (ARC) echelle spectrograph on the 3.5 m telescope at the Apache Point Observatory.

5 Comparison with observations

5.1 Inversion approach

Using the various processes presented in Sect. 2, it is possible to solve the detailed balance equations of C₂ for given (fixed) physical conditions. The free parameters are the temperature, T, of the gas, its number density, n_H, and a scaling factor, G, for the radiation field intensity. If chemical formation excitation is included, then we must also specify the relevant formation rate, k_f.

To solve the “inverse problem”, we define a (χ² like) cost function adapted to excitation diagrams: $$C(T,n_{\mathrm{H}},G,k_f)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{\ln \left(\frac{x_i^{obs}}{g_i}\right)-\ln \left(\frac{x_i^{calc}(T,n_{\mathrm{H}},G,k_f)}{g_i}\right)}{{\tau }_i}\right)}^2,$$(9)

where $x_i^{obs}$ are estimates of the observed fractional populations derived by normalizing⁹ the total observed column density to 1. Here, $x_i^{calc}$ are the fractional populations calculated from the model. For logarithmic variables, the standard maximum likelihood error is ${\tau }_i={\sigma }_i/x_i^{obs}$, where σ_i is the observational error. However, observational fitting procedures give uncertainties that are systematically lower for weak lines than for strong lines. To compensate for this bias, we use ${\tau }_i=\frac{g_i{\sigma }_i}{x_i^{obs}}$ in difficult cases. The number of observed levels is n. It is possible to fix one or several parameters to a prescribed value, and optimize the remaining ones. In Table 1, we report the derived parameters for the three lines of sight that we discuss below and compare with the values derived in Neufeld et al. (2024).

Fig. 8 Observed column densities are scaled to [0 : 1] for comparison with local relative populations xJ. The figure shows a comparison between the present observations and model results for Cyg OB2 12. The observed values are in blue. Model results obtained with the physical parameters derived by Neufeld et al. (2024) and the present study (displayed in Table 1) are reported in green and red, respectively.

5.2 Cyg OB2-12

This line of sight exhibits measurable absorption lines up to J = 32, a very significant improvement over the previously reported data by Neufeld et al. (2024) obtained with ARC echelle data where the highest J is equal to 14, and the data from Hamano et al. (2019) where the highest level is J = 28. However, the J = 30 level has not been detected due to a blend of Q(30) with a diffuse interstellar bands (DIB) feature.

In Fig. 8, we display the excitation diagram obtained from the present observations and the 0D model results derived from the present fit and those reported in Neufeld et al. (2024). We derive a much smaller density, associated with a smaller radiation field scaling. The two models display very similar results up to J = 14 despite these significant differences. We do not compare our results with those derived by Hamano et al. (2019) since their diagnostic was obtained with previous collisional rates.

Standard χ² theory does not allow derivations of accurate error bars here. However, we find that the temperature and chemical rate are well defined, with typical uncertainties of ±2K and ±2 10⁻¹² s⁻¹. The degeneracy between n_H and G is much harder to lift. Fig. 9 shows the cost function computed with the values of T and k_f restricted to their best value. We observe a deep valley, which is not straight, but where n_H/G ≃ 89, and which remains shallow over a large range of densities. The best value reported here could easily shift along that valley under the influence of other physical constraints. We note that using G = 1 would lead to a value of ~89 cm⁻³ for the density, which is close to the value 82 cm⁻³ derived by Neufeld et al. (2024). Here we face the limit of such a 0D model that assumes constant density and temperature. Using a more sophisticated PDR model is probably the best way to incorporate all the necessary ingredients.

Fig. 9 Cost function C for Cyg OB2-12 for T = 34.7 K and kf = 5.6 × 10−12 s−1.

Fig. 10 HD 29647 best fit to (T, nH, G) for a range of kf. Values of kf are color coded. The value selected in Table 1 is shown by a point on the projected curves.

5.3 HD 29647

This line of sight was also reported by Neufeld et al. (2024), where values up to J = 12 are available and considered. Our present observations include levels up to J = 24 but the levels J = 20 and J = 22 are missing, as shown on Fig. 11. Direct minimization of a cost function proved difficult, being very sensitive to computation choices and prone to various degeneracies. So, we fixed parameter k_f to a set of values (from 8 × 10⁻¹³ to 1.5 × 10⁻¹¹ s⁻¹) and for each k_f we derived the optimal set of (T, n_H, G). The results are shown on Fig. 10. For each k_f, we find a well-defined optimum. However, the cost function variations are well below the uncertainties and do not allow us to pick a specific value. We find that T falls within a very narrow range around T = 12 K, and that, also here, there is a well-defined relation between n_H and G, with n_H ≃ 145 × G (cm⁻³).

The set of (degenerate) fits are compared to our observations in Fig. 11. Clearly, an independent argument is needed to lift the degeneracy. Nevertheless, the “best fit” with a negligible chemical excitation is clearly not acceptable. The values shown in Table 1 correspond to the (rather arbitrary) choice k_f = 6× 10⁻¹²s⁻¹.

Fig. 11 Excitation diagram for 13 values of kf going from 0.8 × 10 12 to 15 × 10−12 s−1, using optimal parameters from Fig. 10 (red). The green line corresponds to a negligible chemical excitation.

Fig. 12 Cost function C for T = 38 K and kf = 0 for HD 24534. For all values of nH/G ~ 60, the populations are indistinguishable.

5.4 HD 24534

This is the fiducial case studied by Neufeld et al. (2024) as reported by Sonnentrucker et al. (2007). Since observational column densities are only available up to J = 16, we set k_f = 0. This implies that the degeneracy between n_H and G cannot be lifted, as seen in Fig. 12. The scaling here is n_H ≃ 60 × Gcm⁻³. Fig. 13 shows the proposed fit. One can see that levels J = 0 and J = 16 seem to deviate significantly from the trend of the other levels.

The values reported in Table 1 are only indicative. The associated cloud has an E_B-V ≃ 0.62 (Snow et al. 1998) and is located in the Perseus cloud. We note that C₂ is exposed to a varying radiation field along the line of sight.

6 Summary

In this work we present a detailed model of the excitation of interstellar C₂. Compared to previous work, we extend available molecular data using controlled extrapolations to account for the highest rotational levels of X that are now seen in observations. This includes:

• Extrapolating and deriving rotational electric quadrupole Einstein emission coefficients within the ground electronic state up to J = 34;

• Deriving vibrational 1-0 electric quadrupole Einstein emission coefficients up to v = 1, J = 14;

• Extrapolating and deriving collisional de-excitation rate coefficients with H₂ and He up to J = 34 for ΔJ = 2, 4, 6, and 8.

We have accounted for the electronic radiative transitions between the first eight electronic states of C₂ by using the latest data available in the ExoMol database. These include in particular the A - X Phillips and the a - X intercombination systems. We have also introduced the pumping through the D - X Mulliken system in the UV spectral range by using the data reported in Schmidt & Bacskay (2007) and Krechkivska et al. (2018). We find that their impact is usually negligible.

This all together leads to solving a system of 120 coupled equations including X rotational levels up to J = 34 and levels belonging to a that are located between the rotationally excited levels belonging to X . The contribution of radiative pumping and subsequent radiative decay involving the other levels is computed with the radiative cascade formalism. Our results, limited to these processes, are qualitatively similar to those of Neufeld et al. (2024).

In addition, our coupled equations include the possibility of direct excitation through exothermic chemical formation processes. We show that this excitation mechanism allows us to populate high-lying rotational levels up to the observed values, whereas the standard balance between collisional and radiative processes leads to an unavoidable decline of the excitation diagram for high J values. Our prescription - population proportional to the statistical weight - makes the minimal amount of assumptions, while remaining plausible. It is clear that more studies of this specific process are needed to better quantify the effect.

We show that the lowest rotational levels are mainly sensitive to collisions, which constrains the temperature. Intermediate levels are sensitive to the ratio n_H/G, but not to either of them individually. The highest levels are sensitive to chemical excitation. Their new detection shows that this additional excitation mechanism is relevant for understanding C₂ excitation. If enough J levels are observed, then this allows us to determine a unique set of physical parameters, within the approximation of a 0D model. However, if only the lowest levels are seen, typically up to J = 16, then the degeneracy between density and radiation cannot be lifted. Comparison with previous studies is thus delicate, as most observations are restricted up to J = 16 or less. The availability of determining the population of high J levels then reveals new possibilities. The present trend of our results points to low values of the density, as already found by Neufeld et al. (2024) and resulting from significantly different collisional excitation rate coefficients of C₂ by para and ortho-H₂, compared to previous estimates. We also point out the sensitivity of our results on the kinetic temperature that determines the ortho-to-para ratio of H₂, which further impacts the effect of the collisional excitation of C₂.

This study is a significant step toward a better understanding of C₂ excitation but does not account for the variation of density and temperature in the line of sight driven by the absorption of the radiation field along the cloud. Deriving definitive conclusions on the actual value of the density and cosmic ionization rate may then be premature. Introducing the present excitation model into a full PDR model is our next goal and will allow us to include the constraints from other molecular species as well. Finally, our 0D model also gives the rotational populations in the a state. For the lines of sight that we examined here, their populations are too low to allow for detection.

Fig. 13 HD 24534. Observed column densities are scaled to [0 : 1] for comparison with local relative populations xJ. Model 0D corresponds to the parameters derived in the present study and points in green are obtained with the parameters derived by Neufeld et al. (2024) as shown in Table 1.

Acknowledgements

This work was supported by the Thematic Action “Physique et Chimie du Milieu Interstellaire” (PCMI) of INSU Programme National “Astro”, with contributions from CNRS Physique & CNRS Chimie, CEA, and CNES. We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high performance computing, visualization, and storage resources that have contributed to the results reported within this paper. These results are based on observations obtained with the Habitable-zone Planet Finder Spectrograph on the HET. The HPF team acknowledges support from NSF grants AST-1006676, AST-1126413, AST-1310885, AST-1517592, AST-1310875, ATI 2009889, ATI-2009982, AST-2108512, and the NASA Astrobiology Institute (NNA09DA76A) in the pursuit of precision radial velocities in the NIR. The HPF team also acknowledges support from the Heising-Simons Foundation via grant 2017-0494.

Appendix A C₂ data

A.1 Energy levels

Energies, in cm⁻¹, are in Table A.1.

A.2 Electric quadrupole Einstein coefficients

Einstein emission coefficients are in s⁻¹, and corresponding wavenumbers in cm⁻¹. See Table A.2 and A.3.

A.3 Collision rates extrapolations

Collision rate coefficients with H₂, in cm³ s⁻¹, are in Table A.4 and A.5. Rate coefficients with He are in Table A.6.

Appendix B Cascade coefficients

A cascade coefficient C_nj is the probability that a species on level n > M ends on level j  M after a cascade through any possible path, where M is the number of levels that are involved in the balance equations (1). That probability can be computed using the usual properties of probabilities: it is the sum over all possible intermediate levels of going from n to some level m > j by any path, followed by a direct single transition from m to j. That is: $${\mathcal{C}}_{nj}=\sum _{m=M+1}^n{\mathcal{A}}_{n,m}\frac{A_{m,j}}{A_m}=\sum _{m=M+1}^n{\mathcal{A}}_{n,m}{P}_{m,j}$$(B.1)

where A_m,j is the standard Einstein coefficient and A_m is the inverse radiative lifetime of level m, and we have written $P_{m,j}={\displaystyle \frac{A_{m,j}}{A_m}}$ is the probability to go from n to m by any path and we have ${\mathcal{A}}_{n,n}=1$. These coefficients follow a very similar relation: $${\mathcal{A}}_{n,m}=\sum _{l=m+1}^n{\mathcal{A}}_{n,l}\frac{A_{l,m}}{A_l}=\sum _{l=m+1}^n{\mathcal{A}}_{n,l}{P}_{l,m}$$(B.2)

The only difference is that summation starts at M + 1 for C_{n j} and starts at m + 1 for A_n,m. These coefficients may be computed by recurrence. From any n > M starting from N :

• Set ${\mathcal{A}}_{n,n}=1$

• Then ${\mathcal{A}}_{n,n-1}={\mathcal{A}}_{n,n}{\displaystyle P_{n,n-1}},{\mathcal{A}}_{n,n-2}={\mathcal{A}}_{n,n}{\displaystyle P_{n,n-1}}+{\mathcal{A}}_{n,n-1}{\displaystyle P_{n,n-2}}$ and apply Eq B.2 down to m = M + 1

• Once all pairs (n, m) are known, use Eq B.1 for the final coefficients C_nj.

Appendix C Parameters variations at k_f = 0

To illustrate that only chemical excitation is efficient enough to reproduce the observed population of high J levels, we show here the impact of varying the 3 other parameters (n_H , T and G) around a base model. We use as a reference T = 40 K, n_H = 100 cm⁻³ and G = 1, and we vary each parameter by about 20% around these values. Figs. C.1, C.2 and C.3 show the impact. It is clear that no choice of these parameters is able to give an almost constant value of N_J/g_J for high J, as observed. Chemical excitation succeeds because the state specific formation rate is assumed to be proportional to the statistical weight g_J = 2J + 1. As deexcitation is roughly the same for all high lying levels, this results in a similar value for N_J/g_J.

Table A.4 and A.5 give the rate coefficients shown on Fig. 2. Numbers in parenthesis are powers of 10.

Fig. C.1 Variation of T

Fig. C.2 Variation of nH

Fig. C.3 Variation of G

References

All Tables

All Figures

	Fig. 1 Top: C2 electronic systems included in this computation. “I.C.” stands for “Intercombination” (inspired by Yurchenko et al. 2018). Bottom: Ro-vibrational levels included in the full detailed balance computation. Numerical values are given in Table A.1.
In the text

	Fig. 2 Extrapolation of collisional rate coefficients from Najar et al. (2008) at 30K. Extrapolated rate coefficients are on the right side of the black vertical line. Top: Para-H2, Bottom: Ortho-H2.
In the text

	Fig. 3 Electric quadrupole radiative rotational transition probabilities from van Dishoeck & Black (1982). The positions of J = 22 and 24 are shown by straight lines. The points above are extrapolations.
In the text

	Fig. 4 Adopted ISRF for G = 1 in the range of possible transitions of C2. Colored points displayed on the horizontal line at the bottom of the figure indicate the positions of included transitions coded with the value of Blu Jul from red (strongest) to light green (weakest).
In the text

	Fig. 5 C2 excitation diagram for T = 30K and nH = 50 cm−3 without radiative cascades. Red: Only levels from X included (labeled “X only”); Blue: X and a electronic states, including intercombination transitions (labeled “X + a”).
In the text

	Fig. 6 Effect of nonthermal excitation mechanisms on the weighted fractional populations (see text for the explanation of the labels).
In the text

	Fig. 7 Excitation diagram for a ratio nH/G = 50. All results collapse on a single curve up to J = 16. Note that the density range is much larger than in the cases of Appendix C.
In the text

Fig. 8 Observed column densities are scaled to [0 : 1] for comparison with local relative populations xJ. The figure shows a comparison between the present observations and model results for Cyg OB2 12. The observed values are in blue. Model results obtained with the physical parameters derived by Neufeld et al. (2024) and the present study (displayed in Table 1) are reported in green and red, respectively.

In the text

	Fig. 9 Cost function C for Cyg OB2-12 for T = 34.7 K and kf = 5.6 × 10−12 s−1.
In the text

	Fig. 10 HD 29647 best fit to (T, nH, G) for a range of kf. Values of kf are color coded. The value selected in Table 1 is shown by a point on the projected curves.
In the text

	Fig. 11 Excitation diagram for 13 values of kf going from 0.8 × 10 12 to 15 × 10−12 s−1, using optimal parameters from Fig. 10 (red). The green line corresponds to a negligible chemical excitation.
In the text

	Fig. 12 Cost function C for T = 38 K and kf = 0 for HD 24534. For all values of nH/G ~ 60, the populations are indistinguishable.
In the text

	Fig. 13 HD 24534. Observed column densities are scaled to [0 : 1] for comparison with local relative populations xJ. Model 0D corresponds to the parameters derived in the present study and points in green are obtained with the parameters derived by Neufeld et al. (2024) as shown in Table 1.
In the text

	Fig. C.1 Variation of T
In the text

	Fig. C.2 Variation of nH
In the text

	Fig. C.3 Variation of G
In the text