© The Authors 2026

1 Introduction

The source NGC 7027 is one of the brightest and best-studied young planetary nebulae. It is characterized by a dense, carbon-rich molecular envelope surrounding a hot central star (Sil et al. 2025, hereafter, Paper I and references therein). It serves as a classic example for examining the interactions between intense ultraviolet radiation and molecular chemistry in evolved stellar environments. Despite the harsh ultraviolet radiation field, a rich molecular chemistry is sustained in its outer layers, where vibrationally excited H₂ and CH⁺ have been observed in the photodissociation region (PDR) of NGC 7027. Methylidynium (CH⁺) was initially detected in the planetary nebula NGC 7027 more than two decades ago through far-infrared observations of its pure rotational transitions (Cernicharo et al. 1997). Neufeld et al. (2020) detected nine ro-vibrational emission lines of CH⁺, v = 1 → 0 R(0) − R(3) and P(1) − P(5) for the first time in the same source. Subsequently, Neufeld et al. (2021) reported the detection of four additional ro-vibrational emission lines of CH⁺ v = 1 → 0, P(7) − P(10) transitions along with nine infrared transitions of H₂, comprising the S(8), S(9), S(13), and S(15) pure rotational lines, v = 1 → 0 O(4) − O(7) lines, and the v = 2 → 1 O(5) line. The presence of CH⁺ in this environment is particularly significant because its efficient formation is attributed to the endothermic reaction C⁺(²P) + H₂(¹Σ⁺) → CH⁺ (X¹ Σ⁺) + H(²S), which becomes exothermic when H₂ is rotationally (J ≥ 8, where J is the rotational quantum number) or vibrationally (v ≥ 1, where v is the vibrational quantum number) excited, highlighting the key role of internally excited H₂ in the chemistry within the nebula (Godard & Cernicharo 2013; Faure et al. 2017; Neufeld et al. 2021). If CH⁺ is destroyed before reaching a radiative or collisional equilibrium, the level populations will (partly) reflect its formation conditions. This effect is known as chemical formation pumping and was originally investigated in detail by Godard & Cernicharo (2013). Later, using the first accurate set of state-resolved rate coefficients for CH⁺ rotational excitation (with J ≤ 13) (Faure et al. 2017) and then a larger set with approximate data for CH⁺ ro-vibrational excitation (Neufeld et al. 2021), it was confirmed that chemical pumping indeed has a substantial effect on the excitation of CH⁺ in PDRs.

First identified in absorption in the diffuse interstellar medium (ISM) by Douglas & Herzberg (1941), CH⁺ has subsequently been observed in a variety of astrophysical environments beyond planetary nebulae, including massive star-forming regions (Falgarone et al. 2010), protoplanetary disks (Thi et al. 2011; Zannese et al. 2025), PDRs, such as the Orion bar (Zannese et al. 2025), and external galaxies (Rangwala et al. 2014). We note, in particular, that through the unprecedented spatial resolution and sensitivity of the James Webb Space Telescope (JWST), CH⁺ and CH₃⁺ ions have been detected in emission in the Orion bar and in d203-506 protoplanetary disk regions, heavily irradiated by the Trapezium star cluster (Berné et al. 2023; Zannese et al. 2025).

We report new calculations for the chemistry and excitation of CH⁺ in NGC 7027 based on an accurate set of state-resolved rate coefficients covering the rotational and ro-vibrational excitation of CH⁺. This was made possible through recent theoretical and experimental works providing state-resolved data for the formation, excitation, and destruction of CH⁺, as presented in Appendices A, B, and C. The nonlocal thermodynamic equilibrium (NLTE) analysis of the observed CH⁺ rotational and ro-vibrational emissions is first performed with the code CLOUDY (v23.01; Chatzikos et al. 2023) and then with the code RADEX (van der Tak et al. 2007). Both codes use our newly constructed CH⁺ collision dataset for inelastic transitions. In a second approach, chemical pumping of CH⁺ is implemented in RADEX, and the observed line fluxes are fitted using Markov chain Monte Carlo (MCMC) sampling (Goodman & Weare 2010) performed on the parameter space, namely the kinetic temperature, the atomic hydrogen and electron densities, and the CH⁺ column density. In Sect. 2 we report the predicted results of our NLTE models and compare them with the observed CH⁺ line fluxes. Finally, in Sect. 3, we briefly discuss the implications of our important findings and summarize our concluding remarks.

2 Modeling and results

2.1 CLOUDY modeling and results

We used the 1D photoionization code CLOUDY to model the chemical structure of CH⁺ and to interpret the observed line flux features in the prototypical PDR of NGC 7027, in a similar way as in Paper I (which focused on HeH⁺). To do this, we approximated the elongated nebula to be spherically symmetric and assumed that the dynamically young nebula was in steady-state equilibrium.

2.1.1 Physical conditions

The physical model parameters are listed in Table 1. Our standard model assumes a homogeneous spherical shell of gas surrounding a hot central star, a white dwarf, emitting as a blackbody. This model is computed under the assumption of constant pressure, as in Paper I. The pressure was adjusted by fixing the initial total hydrogen nuclei density [n(H_tot)] at 3.06 × 10⁴ cm⁻³ to obtain a Strömgren sphere of angular radius at 4.6″(~0.02186 pc = 6.74 × 10¹⁶ cm, corresponding to the geometric mean of the semimajor and semiminor axes of the radio continuum emission; Basart & Daub 1987), considering a source distance of 980 pc (Zijlstra et al. 2008). The resulting geometry is a thick, spherical, and closed shell. The default initial gas-phase elemental abundance (relative to the total number of hydrogen nuclei) set for planetary nebulae in CLOUDY modeling is primarily from Aller & Czyzak (1983); Khromov (1989, see Table 2 of Paper I). We only highlight the modifications made to the modeling of NGC 7027 compared to Paper I, while more detailed modeling is presented in Paper I.

Following Lau et al. (2016), we adjusted the dust properties of three independent components with large grains (LGs), very small grains (VSGs), and polycyclic aromatic hydrocarbons (PAHs), constituting exactly 96.6, 1.9, and 1.5% of the total dust by mass, respectively. With these different components of the grain size distribution, we obtained a dust-to-gas mass ratio (M_d/M_g) of ~1/145. A ratio of extinction per reddening of R_V = A_V/E(B − V) = 7.71 and an extinction-to-gas ratio = A_V/N(H) ~ 1.37 × 10⁻²² mag cm² were obtained from our model. For comparison, typical diffuse ISM values are R_V ≈ 3.1 and A_V/N(H) ≈ 5.3 × 10⁻²² mag cm². The significantly higher R_V and lower A_V/N(H) obtained from our model indicate a flatter extinction curve and reduced extinction efficiency per hydrogen atom relative to the diffuse ISM, consistent with a grain population biased toward larger sizes and the dust mass distribution being dominated by large grains, with a relative deficit of small grains.

As a caveat, we note that the sound travel time (~4186 years) in our isobaric model is longer than the kinematic (~600 years; Masson 1989) and dynamic (~1000–2000 years; Santander-García et al. 2012; Ali et al. 2015; Schönberner et al. 2018; Guerrero et al. 2020) ages of the nebula. The formation and destruction timescales (~10⁴ years) of H₂ are also much longer than the age of the system, which makes the steady-state equilibrium assumption practical, but not fully justified. The turbulent velocity dispersion (σ_turb) was set to 6 km s⁻¹ following Neufeld et al. (2021). Finally, our model incorporates the comprehensive model of the H₂ molecule as described by Shaw et al. (2005). In the following, we present the results of our model.

2.1.2 Chemistry

Given the physical conditions described above, the code computed the abundance of each species by solving the chemical reactions under steady-state assumptions. The equilibrium abundance of CH⁺ is dominated by the reactions noted in Table A.1 (see Appendix A), specifically, the formation reaction between C⁺ and H₂ (reaction R1) and destruction reaction with atomic hydrogen (reaction R4) due to their high densities in the PDR. All the rate coefficients of the reactions noted in Table A.1 were updated using recent calculations, except for reaction R5 (for more details, see Appendix A). The other chemical reactions related to CH⁺ followed the default settings in CLOUDY v23.01.

The top panel of Figure 1 shows the equilibrium temperature and density profiles of H, H⁺, H₂, ${\mathrm{H}}_2^{\ast }$ (vibrationally excited H₂), e⁻, C, C⁺, CH⁺, and HeH⁺ as a function of the distance from the star, starting from the illuminated face (~0.0147 pc, where the starlight strikes the cloud) up to a certain depth (A_V = 1.5 mag) where the temperature drops to 80 K (distance of the outer model boundary ~0.0289 pc). In the middle and bottom panels, we present the equilibrium temperature and abundances of species as functions of A_V and the distance (linear and angular) from the star (starting from the ionization front), respectively. The middle and bottom panels are shown explicitly for comparison with the left panel of Fig. 9 in Neufeld et al. (2021), where the authors used the 1D Meudon PDR code (Le Petit et al. 2006) in a plane-parallel setup to model CH⁺ for an isochoric PDR model with a density n(H_tot) = 3 × 10⁵ cm⁻³. We obtained a CH⁺ peak abundance (always with respect to the total number of hydrogen nuclei, unless otherwise stated) of ~2.6 × 10⁻⁷ along the CH⁺ plateau region (~0.0265–0.0275 pc with temperatures of ~800–1500 K). This is lower by a factor of 5 than that obtained (~1.2 × 10⁻⁶) by Neufeld et al. (2021) around temperatures ~1000–2000 K for an isochoric model (see their Fig. 9). It is important to note that the CH⁺ abundance peak from our model (which appears at an angular distance of 5.6″–5.8″ from the central star) shifts toward the outer boundary of the shell as compared to Neufeld et al. (2021), who obtained it around 4.8″–4.9″ (see their left panel of Fig. 9). Thus, our model suggests that the emission of CH⁺ is likely to originate from regions approximately at and extending beyond a spherical shell radius of 5.5″, consistent with the observational angular sizes for pure rotational lines of CH⁺ (see notes in Table 2). The CH⁺ abundance peak plateau falls rapidly at the H/H₂ transition front (the so-called dissociation front where molecular hydrogen number density reaches 50% of the number density of total hydrogen nuclei) around A_V = 1.13 mag. We note that the H/H₂ front is shifted outward (beyond A_V = 1 mag) compared to the findings of Godard & Cernicharo (2013); Neufeld et al. (2021). We identified this effect as resulting from the modified (enhanced) rate coefficient for the destruction of H₂ by C⁺ and from the dust model.

The HeH⁺ density and abundance profiles are also plotted in order to compare them with the findings from Paper I. As in our previous work, the HeH⁺ peak appears near the ionization front (where half of the H⁺ has recombined) with a peak abundance of 1.01 × 10⁻⁸ and a total column density of 1.09 × 10¹² cm⁻², very similar to the results of Paper I. This illustrates the benefit of using the code CLOUDY, which is capable of describing the physical structure of HII regions (from the illuminated face of the cloud to the ionization front) and PDRs.

2.1.3 NLTE calculations

In CLOUDY, the NLTE radiative transfer is treated in 1D geometry within the escape probability formalism (Ferland et al. 2017, and references therein). We employed the closed spherically symmetric geometry (implicit for planetary nebulae) with the sphere expanding (default) option so that the large velocity gradient (LVG), or Sobolev approximation, was employed.

For the spectroscopy and collisional excitation of H₂, the default CLOUDY dataset was employed (Wrathmall et al. 2007). For CH⁺, spectroscopic data (i.e., level energies and radiative rates) were extracted from the EXOMOL database (Tennyson et al. 2016; Pearce et al. 2024) and restricted to the first electronic state (X¹Σ⁺). This set includes a total of 80 ro-vibrational levels, with the highest level (v = 3, J = 11) lying 9495 cm⁻¹ (13 661 K) above the ground state, just below the first excited electronic state, a³Π, at ~1.2 eV (Amitay et al. 1996; Biglari et al. 2014). For inelastic collisions of CH⁺ with hydrogen atoms and free electrons, we combined several accurate sets of theoretical data with approximations based on propensity rules in order to cover all 80 levels, as described in Appendix C.

In Table 2 we compare the observed line fluxes for CH⁺ and H₂ reported by Neufeld et al. (2021), along with those for HeH⁺ and H, He⁺ recombination lines reported by Neufeld et al. (2020) to the predictions of our isobaric CLOUDY model. The line fluxes in CLOUDY are given in units of erg s⁻¹ cm⁻². We converted them into beam-integrated fluxes by scaling to the observed angular sizes (see notes in Table 2) for direct comparison with the observations.

As a goodness diagnostic, we report the mean logarithmic ratio of model-predicted (F_pre) and observed (F_obs) line fluxes, $$\mu =\frac{1}{N}\sum _{i=1}^N\ln \left(\frac{F_{\mathrm{p}\mathrm{r}\mathrm{e}}}{F_{\mathrm{o}\mathrm{b}\mathrm{s}}}\right),$$(1)

where N is the number of transitions. The quantity exp (μ) indicates the geometric mean of the factor by which the model reproduces the observations.

Overall, good agreement is obtained between our CLOUDY model and observations (see orange squares in Fig. 2). On average, the model reproduces observations to within a factor of exp (μ) = 1.28. To be more specific, our model overestimates the CH⁺ pure rotational lines from J = 1 → 0 to J = 4 → 3 by up to a factor of 3, while it underestimates the lines for J = 5 → 4 and J = 6 → 5 by factors of 1.3 and 1.6, respectively. For the CH⁺ ro-vibrational lines, the agreement is within a factor of exp (μ) = 1.22 on average. A good match is thus obtained between our CLOUDY model and observations, suggesting that our updated rate coefficients are reliable and that the chosen H₂ excitation temperature (1900 K) for reaction R1 is appropriate for NGC 7027. It is instructive to compare our results with the best isobaric Meudon PDR model (including CH⁺ chemical pumping) of Neufeld et al. (2021), even though the reactive and inelastic rate coefficients were less accurate. Thus, with a thermal pressure of P/k_B = 3 × 10⁸ K cm⁻³, their model underestimates the pure rotational lines (except J = 1 → 0) and lower-energy ro-vibrational lines, while it overestimates high-excitation ro-vibrational lines of CH⁺ by up to a factor of 2. Thus, globally, the isobaric CLOUDY and Meudon PDR models show comparable agreement with the observations (see blue and orange symbols in Fig. 2). We note that the constant-pressure prescription in CLOUDY refers to total pressure (gas + turbulent + radiation), whereas in the Meudon PDR code, it corresponds to gas thermal pressure alone. Therefore, the numerical pressure values adopted in different codes are not directly equivalent.

In order to better highlight our results, Fig. D.1 illustrates the intrinsic (local) line emissivities of CH⁺ as a function of distance from the star. The rotational and ro-vibrational CH⁺ emission peaks are seen to arise from two distinct physical regions that mainly differ in temperature (~850 versus 1500 K, respectively). We can also observe that the rotational emission peak shows a broader plateau, while the ro-vibrational peak is more sharply peaked. This result suggests that the pure rotational and ro-vibrational lines of CH⁺ probe separate regions in the PDR of NGC 7027, and this might help us to constrain the kinetic temperature gradient just upstream of the H/H₂ transition front.

A comparison between our CLOUDY model-predicted and observed H₂ line fluxes is listed in Table 2 and depicted in Fig. 3. Our isobaric model is found to reasonably explain the observed fluxes of H₂ lines to within a factor of 5.66 on average. As a comparison, the best isobaric Meudon PDR model by Neufeld et al. (2021) reproduces the H₂ ro-vibrational line fluxes within a factor of 4, except for the v = 0 → 0 S(15) and S(13) lines of H₂, which are overestimated by a factor of ~10 and 20, respectively. To facilitate the understanding of the emitting region of the observed H₂ lines, the intrinsic H₂ line emissivities obtained from our CLOUDY model are also plotted in Fig. D.2. Figs. D.1 and D.2 both show that CH⁺ emits from the same region as H₂, meaning that CH⁺ emission is strongly linked to vibrationally excited H₂, as expected.

It should be noted that the dynamical age of NGC 7027 depends on the structures to which it refers. For the inner region traced by the emissions modeled here, it is 1000–2000 years (Santander-García et al. 2012; Ali et al. 2015; Schönberner et al. 2018; Guerrero et al. 2020), while earlier estimates based on multi-epoch observations of the nebula radio emission yield ~600 years (Masson 1989). In the PDR modeled here, H₂ formation is primarily driven by grain-surface reactions and not by the H⁻ channel, even in the warm molecular layer where CH⁺ and H₂ emissions arise at temperatures of ~1500 K. The discrepancy may be due to the time-dependent thermal and chemical evolution of the warm molecular gas, which cannot be captured by the current steady-state approximation. This might explain why our stationary CLOUDY model overpredicts the H₂ emission.

Table 2 shows that for HeH⁺, H I, He II, and He I, our predictions are fully consistent with Paper I. This confirms that our current isobaric model is able to explain most of the observed line fluxes of CH⁺, HeH⁺, H₂, and H, He⁺ recombination lines simultaneously and within a factor of 2–5, with the exception of some ro-vibrational lines of H₂.

Fig. 1 Top panel: equilibrium temperature and densities of H, H+, H2, ${\mathrm{H}}_2^{\ast }$, e−, C, C+, CH+, and HeH+ obtained for an isobaric (total pressure ~2.0 × 109 K cm−3) CLOUDY model with an initial total hydrogen nuclei density n(Htot) = 3.06 × 104 cm−3 as a function of the distance from the star starting from the illuminated face into the depth of the cloud. Middle panel: equilibrium temperature and abundances of species as a function of visual extinction (AV). Bottom panel: equilibrium temperature and abundances of species as a function of the distance (linear and angular) from the star, starting from the illuminated face into the depth of the cloud.

Fig. 2 Ratios of different model-predicted to observed line fluxes (Fpre/Fobs) of rotational and ro-vibrational lines of CH+ as functions of the line names.

Fig. 3 Ratios of CLOUDY model-predicted to observed line fluxes (Fpre/Fobs) of rotational and ro-vibrational lines of H2 as a function of the upper level energies of transitions.

2.2 RADEX modeling and results

In RADEX, the NLTE radiative transfer is solved for a homogeneous and isothermal medium within the escape probability formalism (van der Tak et al. 2007). As with CLOUDY, the LVG approximation (option expanding sphere) was used to model the CH⁺ emission in NGC 7027.

The RADEX parameters for the physical conditions, namely, kinetic temperature T_kin, atomic hydrogen density n(H), and electron density n(e⁻), were extracted from the CLOUDY solution at the CH⁺ abundance peak. As previously discussed in Sect. 2.1.2 and illustrated in Fig. 1, the CH⁺ abundance peak exhibits a plateau. Fig. D.1 shows that the rotational and ro-vibrational CH⁺ lines originate from two distinct physical regions. As noted in Table 3, the main difference between these two regions is the temperature, which increases from 850 K (at the rotational emission peak) to 1500 K (at the ro-vibrational emission peak). Since RADEX considers radiative transfer through a single-zone isothermal medium, radial 1D effects (e.g., varying physical and chemical conditions) can obviously not be captured. However, since the first excited vibrational state (v = 1) of CH⁺ opens at 3942 K above the ground state, it seems more appropriate to use a kinetic temperature above 1000 K in order to efficiently populate by collisions the observed CH⁺ ro-vibrational states. We therefore selected the warmer physical conditions (T = 1507 K) at the ro-vibrational emission peak, allowing the first vibrational state to be more easily populated through collisions via the high-energy tail of the Maxwellian kinetic energy distribution. The last parameter is the CH⁺ column density per unit velocity interval. The total CH⁺ column density was taken as the CLOUDY value (also reported in Table 3), and a line width (FWHM) of 30 km s⁻¹ was selected as representative of the detected lines of CH⁺ (Cernicharo et al. 1997; Neufeld et al. 2021). The incident radiative continuum at the illuminated face of the cloud, derived from the CLOUDY model (and in close agreement with the best-fit dust model of Lau et al. 2016), was also included in RADEX as a background radiation field. We note, however, that the effect of this radiation field on the level populations is a few percent at most, so that radiative pumping of CH⁺ is minor in NGC 7027 (see discussion in Sect. 2.4). Finally, the same spectroscopic and collisional dataset as used in CLOUDY (see Sect. 2.1.3) was used, thus considering inelastic collisions with hydrogen atoms and electrons.

The predicted line fluxes are reported in Table 4. The agreement is better for the ro-vibrational lines (exp (μ) = 0.87) than for the rotational lines (exp (μ) = 0.50), yielding an overall mean agreement of exp (μ) = 0.74. Significant differences between CLOUDY and RADEX are found for CH⁺ pure rotational lines, by factors of ~2–3, whereas for ro-vibrational lines, the RADEX line fluxes are only ~30% lower than those from CLOUDY (see orange and green symbols in Fig. 2, Tables 2 and 4). The discrepancy observed for rotational lines reflects the extended plateau of rotational emissivities primarily, as illustrated in Fig. D.1. Our CLOUDY model indeed predicts a large gradient of the kinetic temperature across this region (just before the H/H₂ front). The inability of RADEX to simultaneously reproduce the observed CH⁺ rotational and ro-vibrational lines with comparable accuracy thus supports the finding that these lines probe distinct regions with different physical conditions. This is corroborated by the much better average factor (exp (μ)) for ro-vibrational lines compared to rotational lines. Finally, we note that the relatively good agreement between CLOUDY and RADEX for the CH⁺ ro-vibrational line fluxes is similar to that observed in Paper I for HeH⁺, whose abundance and line emissivities are even more strongly peaked.

2.3 RADEX/MCMC modeling and results

In this section, another approach that introduces chemical pumping into the NLTE analysis was used to investigate in depth the CH⁺ line excitation. Indeed, because CH⁺ belongs to the class of reactive ions, that is, species that react fast with the dominant colliders (here electrons and H atoms), the chemical source and sink terms must in principle be included when solving the statistical-equilibrium equations. These two terms describe the so-called chemical (or formation) pumping process by which a reactive ion (i.e., with short chemical lifetime) retains memory of its formation conditions. The general consequence is a hotter distribution of rotational and ro-vibrational levels.

The version of RADEX, as modified by Faure et al. (2017), was employed in order to include the state-resolved formation and destruction rate coefficients of CH⁺ for all 80 CH⁺ levels belonging to the ground electronic state X¹Σ⁺ (see Appendix B).

In contrast to the previous (see Sect. 2.2) NLTE RADEX calculations, where the physical conditions and CH⁺ total column density were taken from the CLOUDY results (at CH⁺ ro-vibrational emission peak), RADEX was combined with MCMC sampling in order to determine without a model prescription the physical parameters that best reproduce all 20 lines of CH⁺ detected in NGC 7027. MCMC is a statistical method used to estimate the parameters of a model by iteratively generating random samples from prior probability distributions. To compute the likelihood of the parameters given the observations, we used the publicly available MCMC Python implementation emcee¹ developed by Foreman-Mackey et al. (2013).

As discussed above, the modified version of RADEX was used with the LVG approximation (option expanding sphere), a line width (FWHM) of 30 km s⁻¹, and the background radiation field extracted from CLOUDY. We employed 32 MCMC walkers, each with 4096 steps, to ensure parameter space exploration and convergence. We set the initial priors (noted in Table 5) based on the result obtained from the CLOUDY model near the abundance peak of CH⁺ (see Table 3). The corresponding corner diagram in Fig. E.1 illustrates the posterior distributions and parameter correlations of the four RADEX free parameters (kinetic temperature, atomic hydrogen density, electron density, and CH+ column density) obtained from the RADEX/MCMC sampling. The first 2048 burn-in steps were discarded, and the posterior distributions were constructed from the subsequent 4096 production steps for each of the 32 walkers.

The posterior distributions are found to be well constrained and approximately Gaussian. In particular, the CH⁺ column density is found to be independent of the kinetic temperature and hydrogen atom density. In contrast, N(CH⁺)is anticorrelated with the electron density, reflecting the fact that it is primarily the product of N(CH⁺)and n(e⁻)which is tightly constrained. In practice, the best model gives T_kin = 1258 K, n(H) = 1.62 × 10⁵ cm⁻³, n(e⁻) = 5.89 × 10³ cm⁻³, and N[CH⁺] = 1.48 × 10¹⁴ cm⁻². The kinetic temperature and hydrogen density compare well with the average values derived from CLOUDY around the CH⁺ plateau (see Table 3). In contrast, the electron density is higher by a factor of 10, while the CH⁺ column density is lower by a factor of ~4. As a result, the electron fraction, defined as x_e = n(e⁻)/n(H_tot) and approximated here by x_e ~ n(e⁻)/n(H), is very high, ~3.64 × 10⁻², and much higher than the expected value of ~10⁻³ (i.e., the C⁺ abundance, since most electrons come from the ionization of atomic C; see Fig. 1). On the other hand, the product of N(CH⁺)by n(e⁻)agrees well with that calculated using the CLOUDY values reported in Table 3, which supports the MCMC modeling. In addition, it should be noted that n(e⁻)/n(H) is an upper limit of the actual electron fraction, although the H₂ molecular fraction (f_H2 = 2n(H₂)/[n(H) + 2n(H₂)]) at the CH⁺ abundance peak, as predicted by CLOUDY, is low and typically between 1 and 10% (see Table 3).

The predicted line fluxes are reported in Table 6. The agreement between the MCMC best-fit and the observed CH⁺ line fluxes is better for ro-vibrational lines (to within a factor of exp (μ) = 0.97 on average) than for rotational lines (exp (μ) = 0.50) (see the red symbols in Fig. 2). This result is reminiscent of what was reported previously. Our CLOUDY results suggest that rotational and ro-vibrational lines arise from markedly distinct regions, differing in their temperature and spatial extension. RADEX calculations are therefore intrinsically limited to reproduce all lines (see Sect. 2.2). The overpredicted electron density is likely a consequence of this limitation, namely, the use of a single-zone rather than a one-dimensional approach. It should be noted that removing the ro-vibrational lines from the MCMC fit results in a much smaller electron fraction of ~3 × 10⁻⁵, a much lower kinetic temperature of ~200 K, and a slightly higher CH⁺ column density of ~4 × 10¹⁴ cm⁻². This also comes at the cost of a large underestimation of the ro-vibrational line fluxes. Alternatively, removing the rotational lines from the MCMC fit results in a lower CH⁺ column density of ~3 × 10¹³ cm⁻², and does not substantially improve the electron fraction (~10⁻²).

Fig. 4 Ro-vibrational population diagram of CH+ (v, J) as predicted by RADEX considering the physical condition at ro-vibrational emission peak from the CLOUDY model noted in Table 3.

2.4 Effect of chemical pumping

In order to assess the role of chemical pumping in the CH⁺ excitation, we compare below two RADEX models considering the physical conditions at the ro-vibrational CH⁺ emission peak from the CLOUDY model (see Table 3): one model considers chemical pumping, and in the other, the chemical source and sink terms are ignored in the statistical-equilibrium equations. The results are illustrated in the population diagram of Fig. 4, where the relative population (divided by the degeneracy) of levels with upper energy lower than 8274 K, that is, up to (v = 2, J = 5), is plotted as a function of the upper-level energy. Except for levels (v = 0, J = 1) and to a lesser extent, (v = 1, J = 0), the population of all levels is significantly affected by chemical pumping. This process is found to increase populations by up to a factor of 10. This can be rationalized by comparing the destruction and excitation rates of CH⁺. For the approximate electron fraction x_e ~ n(e⁻)/n(H) = 1.1 × 10⁻³ (at the CH⁺ ro-vibrational emission peak; see Table 3), we expect hydrogen atoms and electrons to dominate the destruction and excitation of CH⁺, respectively (see Appendices B and C). As indicated by our collisional dataset, the destruction state-resolved rate coefficients for H are (at 1300 K) about (4–9) × 10⁻¹⁰ cm³s⁻¹, while the electron-impact excitation rate coefficients from the ground (v, J) = (0, 0) to the first excited level (0, 1) is 1.2 × 10⁻⁶ cm³ s⁻¹, so that excitation of (0, 1) dominates destruction. On the other hand, for higher levels, electron-impact excitation rate coefficients are lower than ~4 × 10⁻⁷ cm³s⁻¹ and destruction prevails. For these levels, the higher populations (with pumping) thus simply reflect their relatively short chemical lifetime. It should be noted that for levels higher than 8274 K, the chemical pumping effect boosts the level populations even further, by several orders of magnitude.

More generally, chemical pumping of CH⁺ should be a major process in most PDR environments because the (state-resolved) destruction of CH⁺ by H or electrons is faster than or comparable to its excitation by these species, with rate coefficients at T_kin = 1000 K of ~5 × 10⁻¹⁰ cm³ s⁻¹ and ~3 × 10⁻⁷ cm³s⁻¹ for destruction, compared to ~10⁻¹⁰ cm³s⁻¹ and ~5 × 10⁻⁷ cm³s⁻¹ for (de)excitation, respectively. In NGC 7027, electrons are the dominant exciting species of CH⁺ because the ionization fraction is large, but excitation by hydrogen atoms would prevail in other environments where x_e ≲ 10⁻⁴. Moreover, radiative excitation was found to play a minor role in NGC 7027 (see Sect. 2.2). This can be understood by comparing the radiative pumping rate to the electron-impact excitation rate, for instance, for the transition (0, 0) → (1, 1): with the radiative continuum flux from CLOUDY at the corresponding wavelength of 3.615 μm (J_ν ~ 1 Jy/nsr), the infrared pumping rate is ~2 × 10⁻⁹ s⁻¹, which is lower by two orders of magnitude than the electron-impact excitation rate of 4.1 × 10⁻⁷ s⁻¹ (at T_kin = 1000 K and n(e⁻) = 5.55 × 10² cm⁻³). Infrared pumping will thus significantly affect the CH⁺ population in environments with much stronger dust continua or much lower densities.

Finally, the line surface brightness variation with wavelength is plotted in Fig. 5. The two spectra correspond to the populations presented in Fig. 4, that is, with and without pumping, for physical conditions at the ro-vibrational emission peak from the CLOUDY model. Table F.1 presents the surface brightness of the most intense CH⁺ ro-vibrational transitions with and without pumping. Chemical pumping increases the most intense lines (those with a surface brightness above 10⁻⁴ erg cm⁻² s⁻¹ sr⁻¹) by factors of ~2–3. For weaker lines, much larger enhancement factors are observed, for instance, for pure rotational lines within v = 1, which at ~30 μm are predicted with line surface brightnesses up to ~2 × 10⁻⁵ erg cm⁻² s⁻¹ sr⁻¹ with pumping, but are lower by one order of magnitude without pumping. The detection of these transitions, as well as ro-vibrational lines of the v = 2 → 1 band at ~4 μm, would provide an even more constraining diagnostic of the physical conditions. We note in this context that six ro-vibrational lines of the v = 2 → 1 band were recently detected with JWST toward the irradiated protoplanetary disk d203-506 (with line intensities of (1–10) × 10⁻⁶ erg cm⁻² s⁻¹ sr⁻¹), where chemical pumping was investigated using a simple excitation model (Zannese et al. 2025).

Fig. 5 CH+ line surface brightness variation with wavelength as predicted by RADEX considering the physical condition at ro-vibrational emission peak from the CLOUDY model noted in Table 3. The left panel shows this without chemical pumping, and the right panel considers chemical pumping. Upper panel: full spectrum in the wavelength range 0–200 μm. Lower panel: most intense ro-vibrational line spectrum in the wavelength range 3–5.5 μm.

3 Discussion and conclusions

The exposure of dense, molecular gas to the UV radiation from the hot central star of the compact young PN NGC 7027 (T ≃ 198 kK, Latter et al. 2000) and to the X-ray emissions from shocked gas (Kastner et al. 2001; Montez & Kastner 2018), represents a rich laboratory for the interplay of radiative and chemical interactions. The conditions between the Strömgren sphere and the inner boundary of the PDR provide a testbed for predictions of the production and destruction rates of reactive ions such as CO⁺ (Bublitz et al. 2023) and CH⁺ (this work). The latter contributes, through a series of two hydrogen abstraction reactions and thus bypassing the path through the slow radiative association of C⁺ and H₂, to the buildup of CH₃⁺ (Smith 1992), which in turn drives the organic chemistry in such environments. In this respect, the first detection of CH₃⁺ in a planetary nebula (NGC 6302, Bhatt et al. 2025, with JWST/MIRI) comes as no surprise, except for the fact that unlike NGC 7027, NGC 6302 is oxygen rich. Conversely, a product of oxygen chemistry, OH, was recently and for the first time detected in a planetary nebula (Ouyang et al. 2024), which is none other than the carbon-rich NGC 7027.

A deeper understanding of these observational findings requires simulations such as those presented in this study. We used the 1D photoionization code CLOUDY and the isothermal single-zone NLTE code RADEX (without and with chemical pumping) to conduct a comprehensive analysis aimed at determining how the observed CH⁺ emissions from NGC 7027 can help us to evaluate our understanding of the CH⁺ chemistry and excitation, and to place constraints on the physical conditions. The analysis was largely based on the recent availability of accurate sets of state-specific rate coefficients for the formation, excitation, and destruction of CH⁺.

Although the limitations are similar as discussed in Paper I, our isobaric CLOUDY model can reproduce observed CH⁺ pure rotational and ro-vibrational line intensities within a factor of 1.3 on average, with a maximum deviation of a factor of 3. The predictions for HeH⁺ and H, He⁺ recombination lines are also entirely consistent with Paper I. On the other hand, our stationary CLOUDY model overpredicts the observed H₂ emissions by up to a factor of 10, perhaps reflecting the short dynamical age of NGC 7027 and the time-dependent chemistry of H₂. Our 1D CLOUDY model also suggests that the rotational and ro-vibrational CH⁺ lines probe physically distinct regions of NGC 7027, with different kinetic temperatures. Our RADEX model with chemical pumping, considering a uniform medium, is indeed unable to reproduce the rotational and ro-vibrational lines intensities simultaneously. The overpredicted electron density (by a factor of 10) obtained from the RADEX/MCMC best-fit model likely reflects the geometric limitation arising from the single-zone approximation.

Conversely, chemical pumping is found to markedly increase the populations of all CH⁺ levels above (v = 0, J = 1). Our RADEX model with chemical pumping thus suggests that the detection of pure rotational lines within v = 1 and ro-vibrational lines of the v = 2 → 1 band would provide additional constraints on the physical conditions, especially the kinetic temperature gradient, just upstream of the H/H₂ transition region. We conclude that it is highly desirable to enable CLOUDY to consistently couple subsets of the underlying chemical reaction network to the statistical equilibrium calculations. It will also be interesting to implement the current set of state-resolved inelastic and reactive rates in a 1D complete PDR model, such as that of Neufeld et al. (2021).

Data availability

The new extended CH⁺ ro-vibrational inelastic dataset underlying this article will be made available in the Excitation of Molecules and Atoms for Astrophysics (EMAA) database, https://emaa.osug.fr and https://dx.doi.org/10.17178/EMAA (Faure et al. 2025).

Acknowledgements

The authors thank the anonymous referee for the valuable suggestions and constructive feedback, which have significantly improved the quality of this work. M.S. acknowledges financial support from the European Research Council (consolidated grant COLLEXISM, grant agreement ID: 811363). This research was carried out within the framework of IPAG, UMR 5274, with support from a CNRS postdoctoral fellowship. M.S. further acknowledges support from the National Science and Technology Council, Taiwan (Grant Nos. NSTC-114-2811-M-007-017, NSTC-111-2112-M-007-014-MY3, NSTC-113-2639-M-A49-002-ASP, and NSTC-113-2112-M-007-027), and from the Postdoctoral Incentive Award of the Research Talent Resource Center at National Tsing Hua University. T.G.L. acknowledges funding from Grant No. MICIU/AIE/10.13039/501100011033 PID2024-155666NB-I00. J.F. acknowledges support from the USA National Science Foundation award number AST-2303895. This research made use of CLOUDY (Chatzikos et al. 2023), RADEX (van der Tak et al. 2007), emcee (Foreman-Mackey et al. 2013), astropy (http://www.astropy.org; Astropy Collaboration 2022), corner (Foreman-Mackey 2016), matplotlib (Hunter 2007), numpy (Harris et al. 2020), UKRmol+ 3.2 (for the electron-CH⁺ scattering calculations; Mašín et al. 2020), and Psi4 (for the underlying CH⁺ electronic structure calculations; Smith et al. 2020). This document was prepared using the collaborative tool Overleaf available at: https://www.overleaf.com/.

References

Appendix A Updated reaction rates

The CH⁺ chemical network from CLOUDY version c23.01 was updated for five among the six dominant formation and destruction reactions R1–R6, as listed in Table A.1. For the hydrogenation of C⁺(²P) by H₂ (reaction R1), we follow CLOUDY’s distinction between two forms of molecular hydrogen: H₂ represents the ground vibrational state (v = 0) and ${\mathrm{H}}_2^{\ast }$ refers to vibrationally excited states (v ≥ 1). CLOUDY defines ${\mathrm{H}}_2^{\ast }$ as all levels in the X electronic state with an energy greater than 4100 cm⁻¹. In practice, we assumed that H₂ corresponds to the first 7 ro-vibrational pure rotational levels (v = 0, J = 0 → 6), whose energies are significantly below the endothermicity $\mathrm{\Delta }{H}_0^{\mathrm{o}}=4620\mathrm{K}$ of the reaction C⁺(²P) + H₂(¹Σ⁺) → CH⁺ (X¹Σ⁺) + H(²S) (Hierl et al. 1997). Higher levels were attributed to ${\mathrm{H}}_2^{\ast }$, that is, the form of molecular hydrogen which is reactive with C⁺(²P). A uniform excitation temperature of 1900 K was adopted for all levels of molecular hydrogen, which was found to reproduce the rotational diagram for the observed transitions of H₂ fairly well (Neufeld et al. 2021). The state-to-state rate coefficients of González-Lezana et al. (2026) (see Appendix B) were finally employed to derive the rate coefficients for H₂ and ${\mathrm{H}}_2^{\ast }$ listed in Table A.1. For the photoionization of CH (reaction R2), we employed the A_V-dependent rate coefficient from Heays et al. (2017) with two different product channels as recommended by KIDA². For the dissociative recombination (DR) of CH⁺ (reaction R3), the rate coefficient measured by Paul et al. (2022) for the ground state of CH⁺ (v = 0, J = 0) was used. It should be very similar to the thermal rate coefficient given the weak J-specific dependence of the DR process (see Appendix B). For the reaction of CH⁺ with H (reaction R4), we employed the thermal time-independent quantum mechanical (TIQM) rate coefficient computed by Faure et al. (2017), which agrees well with available experimental measurements above 50 K (see Appendix B). For the reaction of CH⁺ with H₂ (reaction R5), the default CLOUDY value was employed, as recommended by KIDA and UMIST³. Finally, for the photodissociation of CH⁺ (reaction R6), the A_V-dependent rate coefficient was taken from Heays et al. (2017). For the reactions R1, R3, and R4, the rate coefficients were fitted using a standard modified-Arrhenius equation. We note that our fits are reliable over the kinetic temperature ranges 100 − 3000 K for reaction R1, 10 − 10 000 K for reaction R3, and 10 − 3000 K for reaction R4, and they should be accurate to within ~20%. We further point out that Paper I contains a typographical error in the thermal rate coefficient for reaction R1 between He⁺ and H; the corrected form is noted here in Table A.1.

Appendix B New CH⁺ ro-vibrational data for reactive collisions

This dataset was employed in the RADEX/MCMC calculations (see Sect. 2.3).

B.1 Hydrogenation of C⁺(²P)

For the reaction of C⁺(²P) with H₂(v, J), which produces CH⁺(v′, J′) + H, we employed the extensive set of rate coefficients computed recently by González-Lezana et al. (2026) using the statistical quantum method (SQM). The use of statistical methods is justified for this system because the intermediate complex CH₂⁺ forms a deep potential well (Konings et al. 2021; González-Lezana et al. 2026). The SQM calculations were performed for the first 31 levels of H₂ (up to (v, J) = (1, 10) at 14 223 K above the ground state) and for the first 42 levels of CH⁺ (up to (v′, J′) = (2, 5) at 8275 K above the ground state). They were compared to results obtained with a different statistical approach (Konings et al. 2021) and to the results of Faure et al. (2017), based on the more accurate but expensive time-dependent wave-packet (TDWP) method, and a good agreement (within a factor of 2) was observed at the state-to-state level. This is illustrated in Fig. B.1 where we have plotted the occupation probability of CH⁺ following the reaction of C⁺(²P) with the excited state of H₂ (v′, J′) = (2, 0) at T_kin = 1000 K⁴. We note that a good agreement was also observed between the SQM calculations and experiments performed with H₂ (v = 0, 1) at different temperatures, see e.g., Fig. 1 of González-Lezana et al. (2026). The full set of SQM rate coefficients includes 1194 state-to-state transitions and covers the kinetic temperature range from 30 to 3000 K. We note that the original reactive rate coefficients of González-Lezana et al. (2026) are limited to the temperature range 30 − 1500 K to ensure full convergence. They were extended here up to 3000 K to better cover the physical conditions of NCG 7027, with limited impact on accuracy (see González-Lezana et al. 2026).

In order to cover all 80 CH⁺ levels in its ground electronic state X¹Σ⁺, we have also extrapolated the data of González-Lezana et al. (2026) (that is, above (v′, J′) = (2, 5)) using a simple formulation similar to that of Neufeld et al. (2021) (see their Appendix B.2), where a formation temperature is introduced additionally. The (normalized) probability p(v′, J′) of forming CH⁺ in a level (v′, J′) is defined as: $$p(v^{\mathrm{\prime }},J^{\mathrm{\prime }})=\{\begin{array}{ll}\frac{(2J+1)\exp (-E/T_{\text{form\hspace{0.17em}}})}{\sum (2J+1)\exp (-E/T_{\text{form\hspace{0.17em}}})}& \text{\hspace{0.17em}if\hspace{0.17em}}E<0\\ \frac{(2J+1)\exp (-E/T_{\text{kin\hspace{0.17em}}})}{\sum (2J+1)\exp (-E/T_{\text{kin\hspace{0.17em}}})}& \text{\hspace{0.17em}if\hspace{0.17em}}E\ge 0,\end{array}$$(B.1)

where $E=\mathrm{\Delta }{H}_0^{\mathrm{o}}+E_{{\mathrm{C}\mathrm{H}}^{+}}({\mathit{v}}^{\mathrm{\prime }},J^{\mathrm{\prime }})-E_{{\mathrm{H}}_2}(\mathit{v},J)$ is the activation energy, including the reaction endothermicity $\mathrm{\Delta }{H}_0^{\mathrm{o}}=4620\mathrm{K}$ (Hierl et al. 1997) and E_CH⁺(v′, J′) and E_H2(v, J) are the ro-vibrational energies of CH⁺ and H₂, respectively. The formation temperature T_form was taken as $2/3\times (E_{{\mathrm{H}}_2}(v,J)-\mathrm{\Delta }{H}_0^{\mathrm{o}})$, as observed and recommended by Faure et al. (2017). We note that in the measurements of the deuterated variants of the H₂ + ${\mathrm{H}}_2^{+}$ reaction, about two-thirds of the reaction exothermicity was also found in the internal ro-vibrational energies of the (triatomic) product ions (Merkt et al. 2022, although their population does not necessarily follow a Boltzmann distribution). As can be seen in Fig. B.1, the above extrapolation scheme (Extrap) reproduces the SQM probabilities to within a factor of 2 for the first 42 levels of CH⁺.

B.2 Dissociative recombination of CH⁺ with electrons

For the DR of CH⁺ (v, J) with electrons, we have computed new cross sections using the theoretical treatment of Forer et al. (2024) combining the R-matrix method with the multichannel quantum defect theory (MQDT), adopting here a partial wave basis l = 0 − 4 and electron energies up to 1 eV. Rate coefficients were obtained for the CH⁺ levels (v = 0, J = 0 − 11) and (v = 1, J = 0 − 11). As in Forer et al. (2024), these rate coefficients could be compared to state-resolved measurements made at the Cryogenic Storage Ring (CSR) by Paul et al. (2022). These authors were able to determine the J-specific rate coefficient for the levels J = 0, 1 of CH⁺ in its ground electronic and vibrational state, and for kinetic temperatures between 10 and 40 000 K. Our new DR rate coefficients for v = 0, J = 0–1 agree with those of Paul et al. (2022) to within a factor of ~2–3 in the range 10–3000 K, as illustrated in Fig. B.2 for J = 0. We can also observe that the rotational effects are weak, with J-specific DR rate coefficients for CH⁺(v = 0) differing by less than a factor of 3 at 1000 K. Moreover, the DR rate coefficients of CH⁺(v = 1) are found to exceed those of CH⁺(v = 0) by a small factor (<3), and rotational effects are even less significant. In practice, the MQDT rate coefficients for CH⁺(v = 0, J = 11) were used for all levels (v = 0, J > 11) and the rate coefficients for (v = 1, J = 11) were employed for all levels (v = 1, J > 11) and (v > 1, J). We thus selected only theoretical DR rate coefficients, rather than the experimental ones for (v = 0, J = 0–1), to ensure consistency among the (v, J)-specific DR rate coefficients and with the excitation rate coefficients described below.

Fig. B.1 Occupation probability of the nascent CH+ product from the reaction of C+(2P) with H2 in the ro-vibrationally excited state (v, J) = (2, 0) as function of the CH+(v′, J′) level energy and at a kinetic temperature Tkin = 1000 K. The SQM calculations of González-Lezana et al. (2026) are compared to the TDWP results of Faure et al. (2017) and to the simple extrapolation scheme (Extrap) defined by Eq. B.1. See text for details. The rotational ground states in the rotational manifold of v′ = 0, 1, 2 are labeled for clarity.

Fig. B.2 State-selected DR rate coefficients of CH+(v, J). The present MQDT calculations are denoted by colored solid lines. The thick black dashed line gives the experimental rate coefficients from Paul et al. (2022) for J = 0.

B.3 The reaction of CH⁺ with H atoms

For the reaction of CH⁺ with H (which produces C⁺(²P) + H₂), we used TIQM calculations of Faure et al. (2017) for CH⁺ in levels (v = 0, J = 13) and the SQM calculations of González-Lezana et al. (2026) for higher levels up to (v = 2, J = 5), consistent with the calculations described above for the reverse reaction between C⁺(²P) and H₂. As for this latter reaction, this set of state-resolved rate coefficients (summed over the final ro-vibrational states of H₂) was found to be in good agreement with the more accurate TIQM calculations of Faure et al. (2017), typically within a factor of 2. The SQM thermal rate coefficient thus matches well with the thermal TIQM value and also with various measurements, except below ~ 50 K where a strong decrease of the rate coefficient was observed in the ion-trap experiment of Plasil et al. (2011). This surprising result was interpreted by these authors as a dramatic loss of reactivity of rotationally cold CH⁺ ions at very low temperatures. As recently discussed by del Mazo-Sevillano et al. (2025), the vibrational excitation of CH⁺ strongly reduces reactivity at low temperatures, so that vibrational thermalization was perhaps not fully achieved in the ion-trap experiment (see Fig. 9 of del Mazo-Sevillano et al. 2025). A similar result was observed by González-Lezana et al. (2026). For levels not included in the SQM set of González-Lezana et al. (2026), that is, for levels above CH⁺(v = 2, J = 5), we used an average temperature-independent rate coefficient of 4 × 10⁻¹⁰ cm³ s⁻¹. We note that rotational and vibrational effects are weak for kinetic temperatures above ~500 K, with ro-vibrational rate coefficients lying in the range (4–5) × 10⁻¹⁰ cm³ s⁻¹. This range is ~ 2 − 3 orders of magnitude lower than the above DR rate coefficients so that hydrogen atoms will dominate the destruction of CH⁺ in regions where n(e)/n(H) < 2 × 10⁻³.

Appendix C New CH⁺ ro-vibrational data for inelastic collisions

This dataset was employed in the CLOUDY and RADEX calculations (see Sect. 2.1.3 and 2.2).

C.1 Electron collisions

For inelastic collisions of CH⁺(v, J) with electrons, new cross sections were computed using the theoretical treatment of Forer et al. (2024), adopting a partial wave basis l = 0 − 4 and electron energies up to 1 eV, as for DR. Pure rotational and ro-vibrational rate coefficients were obtained for the CH⁺ levels (v = 0, J = 0 − 11) and (v = 1, J = 0 − 11) − all calculated as described by Forer et al. (2024). This method includes the lowest two electronically excited states of the ion, which produce low-energy resonances in the scattering phases of the ground electronic state and can significantly increase the vibrational coupling between different vibrational levels. For transitions with ΔJ = ±1, the limited partial-wave basis of the scattering calculations is corrected via Born closure with the Coulomb-Born (CB) approximation, as in Hamilton et al. (2016). The CB approximation accounts for the larger partial waves that are omitted in the R-matrix scattering basis; these are important due to strong dipole-induced coupling between partial waves that satisfy Δl = ±1. The CB approximation is known to be reliable for strongly dipolar ions like CH⁺ (1.68 D), but does not include resonances in the scattering phases and, therefore, is only used to supplement the limited partial-wave basis of the R-matrix calculations. In fact, the CB approximation on its own was found to underestimate the ro-vibrational rate coefficients by three orders of magnitude overall.

As in Forer et al. (2024), our pure rotational rate coefficients could be compared to state-selected measurements made at the CSR by Kálosi et al. (2022). These authors were able to evaluate the rate coefficients for the excitation J = 0 → J = 1 of CH⁺ in its ground electronic and vibrational state, for kinetic temperatures up to 200 K. Our new excitation rate coefficients agree with those of Kálosi et al. (2022) to within 50% across the 10 − 200 K temperature, as illustrated in Fig. C.1, which is consistent with the experimental one-sigma uncertainty. We can also observe that the dominant transitions are those with ΔJ = 1, 2 for rotational excitation and ΔJ= 0, 1, 2 for ro-vibrational excitation. Another important result is that ro-vibrational excitation is dominated by transitions with Δv = 1 (see Fig. 2 of Forer et al. 2024) and that the corresponding rate coefficients do not significantly depend on the initial state v, especially above 100 K. As a result, the pure rotational rate coefficients for de-excitations ΔJ = −1, −2 of CH⁺(v = 0, J = 11) and (v = 1, J = 11) were used for all levels (v = 0, J > 11), (v = 1, J > 11) and (v ≥ 1, J), respectively. Similarly, the ro-vibrational rate coefficients for (de-)excitations Δv = −1; ΔJ = 0, ±1, ±2 of CH⁺(v = 1, J = 9) were used for all levels (v = 1, J > 11) and (v > 1, J). The full set of CH⁺ + e⁻ inelastic rate coefficients is available for temperatures between 10 and 3000 K.

Fig. C.1 State-to-state excitation rate coefficients from the ground state of CH+(v = 0, J = 0) to the lowest three and four rotational and ro-vibrational excited states (v′, J′), respectively. The present MQDT calculations (CB corrected) are denoted by colored solid lines. The thick black dashed line gives the experimental rate coefficients from Kálosi et al. (2022) for v′ = 0, J′ = 1.

C.2 Collision with hydrogen atoms

For inelastic collisions of CH⁺(v, J) with hydrogen atoms, we used the TIQM calculations of Faure et al. (2017) for CH⁺(v = 0, J = 0 − 13) and the SQM rate coefficients computed by González-Lezana et al. (2026) for levels above (v = 2, J = 5), consistent with the calculations used for the CH⁺ + H reactive channel described above. We note that since CH⁺ + H has a statistical behavior (González-Lezana et al. 2026, and references therein), inelastic collisional propensity rules are different from non-statistical systems like, e.g., HeH⁺ (Paper I, and references therein). In particular, for a specific initial ro-vibrational state (v, J), the ro-vibrational rate coefficients for a transition (v, J) → (v′, J′) do not depend on v′ but only on J′. Thus, for levels higher than (v = 2, J = 5), we extrapolated the data of González-Lezana et al. (2026) by applying a temperature-independent de-excitation rate coefficients of 5 × 10⁻¹¹ cm³ s⁻¹ for all rotational and ro-vibrational transitions with |ΔJ| ≤ 5. This crude prescription should be accurate to within a factor of 5 − 10. We note that the rate coefficients for the rotational (ro-vibrational) excitation of CH⁺ due to H collisions are typically four (three) orders of magnitude smaller than those due to electron collisions, meaning that hydrogen atoms will dominate the CH⁺ excitation if n(e)/n(H) < 10⁻⁴. The full set of CH⁺ + H inelastic rate coefficients is available for temperatures between 30 and 3000 K. We note that the original inelastic rate coefficients of González-Lezana et al. (2026) are limited to the temperature range 30 − 1500 K to ensure full convergence. They were extended here up to 3000 K to better cover the physical conditions of NCG 7027, with limited impact on accuracy (see González-Lezana et al. 2026).

Appendix D Emissivity profiles of CH⁺ and H₂

Figures D.1 and D.2 display the intrinsic line emissivities of CH⁺ and H₂, respectively, as a function of the distance from the star.

Fig. D.1 Emissivity profile of CH+ pure rotational lines and ro-vibrational lines as a function of distance.

Fig. D.2 Emissivity profile of the observed H2 lines as a function of distance.

Appendix E RADEX/MCMC fitting results

The corner diagram resulting from the RADEX/MCMC analysis of CH⁺ line fluxes is presented in Fig. E.1.

Fig. E.1 Corner plot obtained from the RADEX/MCMC exploration for CH+, taking into account all detected rotational and ro-vibrational transitions. The diagonal panels show the marginalized posterior probability distributions of Tkin (linear scale), n(e−), n(H), and N(CH+) (log10 scale), while the off-diagonal panels display the corresponding pair-wise correlations. The orange lines and markers indicate the maximum posterior estimate, and the vertical dashed lines in each diagonal panel mark the 16th and 84th percentiles (i.e., the 1σ credible intervals). The contours in the 2D panels enclose 39.3% and 86.4% of the posterior probability, corresponding to the 1σ and 2σ confidence regions for a two-dimensional Gaussian distribution.

Appendix F CH⁺ line surface brightness prediction

The surface brightness values of the strongest CH⁺ ro-vibrational transitions in the presence and absence of chemical pumping are noted in Table F.1. The values are predicted by RADEX considering the physical condition at the ro-vibrational emission peak and CH⁺ total column density obtained from the CLOUDY model noted in Table 3. A line width (FWHM) of 30 km s⁻¹ and the incident radiative continuum at the illuminated face of the cloud, obtained from the CLOUDY model, were used as the background radiation field.

All Tables

All Figures

Fig. 1 Top panel: equilibrium temperature and densities of H, H+, H2, ${\mathrm{H}}_2^{\ast }$, e−, C, C+, CH+, and HeH+ obtained for an isobaric (total pressure ~2.0 × 109 K cm−3) CLOUDY model with an initial total hydrogen nuclei density n(Htot) = 3.06 × 104 cm−3 as a function of the distance from the star starting from the illuminated face into the depth of the cloud. Middle panel: equilibrium temperature and abundances of species as a function of visual extinction (AV). Bottom panel: equilibrium temperature and abundances of species as a function of the distance (linear and angular) from the star, starting from the illuminated face into the depth of the cloud.

In the text

	Fig. 2 Ratios of different model-predicted to observed line fluxes (Fpre/Fobs) of rotational and ro-vibrational lines of CH+ as functions of the line names.
In the text

	Fig. 3 Ratios of CLOUDY model-predicted to observed line fluxes (Fpre/Fobs) of rotational and ro-vibrational lines of H2 as a function of the upper level energies of transitions.
In the text

	Fig. 4 Ro-vibrational population diagram of CH+ (v, J) as predicted by RADEX considering the physical condition at ro-vibrational emission peak from the CLOUDY model noted in Table 3.
In the text

Fig. 5 CH+ line surface brightness variation with wavelength as predicted by RADEX considering the physical condition at ro-vibrational emission peak from the CLOUDY model noted in Table 3. The left panel shows this without chemical pumping, and the right panel considers chemical pumping. Upper panel: full spectrum in the wavelength range 0–200 μm. Lower panel: most intense ro-vibrational line spectrum in the wavelength range 3–5.5 μm.

In the text

Fig. B.1 Occupation probability of the nascent CH+ product from the reaction of C+(2P) with H2 in the ro-vibrationally excited state (v, J) = (2, 0) as function of the CH+(v′, J′) level energy and at a kinetic temperature Tkin = 1000 K. The SQM calculations of González-Lezana et al. (2026) are compared to the TDWP results of Faure et al. (2017) and to the simple extrapolation scheme (Extrap) defined by Eq. B.1. See text for details. The rotational ground states in the rotational manifold of v′ = 0, 1, 2 are labeled for clarity.

In the text

	Fig. B.2 State-selected DR rate coefficients of CH+(v, J). The present MQDT calculations are denoted by colored solid lines. The thick black dashed line gives the experimental rate coefficients from Paul et al. (2022) for J = 0.
In the text

	Fig. C.1 State-to-state excitation rate coefficients from the ground state of CH+(v = 0, J = 0) to the lowest three and four rotational and ro-vibrational excited states (v′, J′), respectively. The present MQDT calculations (CB corrected) are denoted by colored solid lines. The thick black dashed line gives the experimental rate coefficients from Kálosi et al. (2022) for v′ = 0, J′ = 1.
In the text

	Fig. D.1 Emissivity profile of CH+ pure rotational lines and ro-vibrational lines as a function of distance.
In the text

	Fig. D.2 Emissivity profile of the observed H2 lines as a function of distance.
In the text

Fig. E.1 Corner plot obtained from the RADEX/MCMC exploration for CH+, taking into account all detected rotational and ro-vibrational transitions. The diagonal panels show the marginalized posterior probability distributions of Tkin (linear scale), n(e−), n(H), and N(CH+) (log10 scale), while the off-diagonal panels display the corresponding pair-wise correlations. The orange lines and markers indicate the maximum posterior estimate, and the vertical dashed lines in each diagonal panel mark the 16th and 84th percentiles (i.e., the 1σ credible intervals). The contours in the 2D panels enclose 39.3% and 86.4% of the posterior probability, corresponding to the 1σ and 2σ confidence regions for a two-dimensional Gaussian distribution.

In the text