Multistate coupled diabatic neural network potential for the quantum non-adiabatic photofragmentation of CH₂⁺

Abstract

Tracking the complex non-adiabatic transitions in far-ultraviolet photodissociation demands highly accurate diabatic potential energy matrices (PEMs) across numerous excited states. To address this, we introduce a fully automated diabatization method that leverages artificial neural networks to fit PEMs. Our approach divides the PEM into a physically grounded zeroth-order diagonal term, which is then corrected by a neural network matrix to capture electronic couplings. By enforcing symmetry constraints on off-diagonal elements and sharing degenerate diabatic states between the A′ and A″ irreducible representations, the diabatization process becomes completely automatic. We validate this method using time-dependent wavepacket calculations to simulate the photodissociation of CH₂⁺, incorporating relevant states up to ≈13.6 eV. Finally, we compute partial cross-sections for all fragmentation channels—including total and partial fragmentation yielding CH⁺, CH, H₂, and H₂⁺ diatoms—revealing a notably high cross-section for the formation of the CH radical.

---

1 Introduction

Photodissociation plays a fundamental role in outer atmospheres, on Earth and exo-planets, as well as in the borders of molecular clouds and protoplanetary disks,^1–4 by controlling the abundance of small molecules. Those environments are irradiated by high-energy ultraviolet photons producing a large variety of radicals and ions, which in turn trigger the formation of more complex molecules. Hence, it is not only crucial to know the total photodissociation rates, but also the branching ratios for the formation of photo-products in different electronic states.

The fragmentation dynamics at high energy is governed by non-adiabatic transitions among many excited electronic states. When these states cross, leading to conical intersections (CI), the dynamical simulations become complicated, since the non adiabatic couplings (NACs) present singularities.^5,6 This problem is frequently overcome by a unitary transformation to a diabatic electronic representation^7,8 where the kinetic adiabatic couplings vanish. However, this is in general not possible for all degrees of freedom,⁹ and the term quasi-diabatization is commonly used. Several diabatization procedures exist and they are classified as derivate-, property- or energy-based methods, depending on the information they use.¹⁰ One common approach is the regularization,¹¹ i.e. the elimination of the singularities appearing at conical intersections.

Block diagonalization^12,13 and effective Hamiltonian¹⁴ methods use reference geometries to generate a local description of diabatic Hamiltonians and are successfully applied to study photodissociation spectra. H₃⁺,¹⁵ C₆H₆⁺ (ref. 16) and C₄N₂H₄¹⁷ are representative examples of diabatization, using derivative-based methods, of large systems including four or five electronic states, using normal coordinates thus allowing the accurate description of photoabsorption spectra using quantum time dependent methods, but not fully follow the dissociation dynamics. However, in many cases there are multiple crossings, not only in the Frank–Condon region but also close to products in different rearrangement channels, that need to be considered to properly describe the branching ratios in photodissociation. The problem becomes more complex when dealing with large energy intervals, as in astrochemistry, because a large number of electronic states are needed to describe the photofragmentation rate in all the UV radiation field.

The application of machine learning techniques to produce diabatic potential energy matrices (PEM) is being found to be extremely useful as has been already proven in different systems.^18–21 For two-state systems fully automatic diabatization methods (meaning no diabatic information is used) can provide reasonable PEM from pure adiabatic information,^22–24 leveraging the symmetry restrictions in the couplings of the PEM. As the number of states grows larger semi-automatic diabatization methods, as the family of diabatization by deep neural networks (DDNN²⁵), exists, where minimal diabatic information is required.^26,27 Other approaches use machine learning methods with clustering and regression techniques to generate smooth and well ordered diabatic states with property-based diabatization.²⁸ In this work we present a fully automatic diabatization method applicable to many electronic states over a large energy interval for triatomic systems. For that we first construct an initial PEM model for describing the products in different rearrangements, which are correlated by symmetry considerations allowing to impose an ordering of the diabatic states. A many-body term expressed as an artificial neural network is added to the zeroth order model to correct the short distance region and include the couplings among the diabatic states with symmetry restrictions.

The diabatization method is applied to the study of the photodissociation of CH₂⁺ cation, whose ground state presents important Renner–Teller effects due to its ²Π character²⁹ and it has multiple ionic crossings in the CH⁺/CH and H₂/H₂⁺ products rearrangement channels, shown in Fig. 1, thus providing a good benchmark. This system is a missing stone in the hydrogenation chain CH_n⁺ + H₂ → CH_n+1⁺ + H, giving rise to CH₃⁺ which is considered the origin of hydrocarbons in the interstellar medium (ISM).^30,31 CH⁺ is one of the first molecules observed in the ISM,³² and since then has been observed in a great variety of interstellar and circumstellar environments. On the other hand, CH₃⁺ was observed very recently in a protoplanetary disk with the James Webb Space Telescope.³³ Thus the missing detection of CH₂⁺ is fundamental to corroborate the above mentioned hydrogenation chain as the formation mechanism of CH₃⁺. Formation and destruction rates are needed to properly model the abundance of these ions. The photodissociation cross sections have been calculated for CH⁺ (ref. 34 and 35) and CH₃⁺ (ref. 36) using quantum dynamical methods, while for CH₂⁺ only a vertical excitation model from the ground equilibrium geometry is available.³⁷

Fig. 1 Top panel: Diatomic diabatic potential energy curves for the [CH + H]⁺ asymptote. Bottom panel: Diatomic diabatic potential energy curves for the [C + H₂]⁺ asymptote. The zero of energy is placed at the minimum of the CH₂⁺ in the ground electronic state X²A′.

The final goal of this work is to study the non-adiabatic photofragmentation dynamics of CH₂⁺ in a new diabatic model composed of 16 coupled electronic states, and it is organized as follows. First, the new diabatization method is presented and applied to CH₂⁺. Next, the results on the photodissociation dynamics obtained with a quantum wave packet treatment are discussed. Finally, some conclusions are extracted.

2 Diabatization

Machine learning diabatization methods often rely on some diabatic information to impose the correct order of the different states. In this work we present a fully automatic diabatization method based on a many body expansion of the potential energy. It is applied to produce a diabatic potential energy matrix (PEM) of 16 states of the CH₂⁺ system, 8 in each A′ and A″ irreducible representations of the C_s point group. With this number of states we properly include the Frank–Condon region of the CH₂⁺ and possible dissociation channels up to ≈13.6 eV. While the diabatization method is automatic, it still requires a very detailed analysis of the electronic structure of the system.

2.1 Diatomic fragments

The diatomic fragments in the [CH + H]⁺ and [C + H₂]⁺ channels in Fig. 1 present a great variety of electronic states, bound and dissociative, in which the diatomic fragment can be neutral or ionized. When performing the calculations for the whole triatomic system in the A′ or A″ representations the adiabatic states present many crossings, related for instance with charge transfers between the distant partners. These adiabatic states are first diabatized by performing independent calculations of each diatomic fragment, with different charges and in different irreducible representations. Fig. 2 displays the diabatized diatomic curves for both CH/CH⁺ and H₂/H₂⁺ systems.

Fig. 2 Representation of the diabatic diatomic curves included in the U⁰ model and correlation with the total dissociation asymptotes.

2.2 Correlation diagram

The different electronic states in all the rearrangement channels are connected with the lowest energy total fragmentation asymptotes: C⁺(²P) + H(²S) + H(²S), C(³P) + H(²S) + H⁺ and C(¹D) + H(²S) + H⁺. A full description of the 22 states that correlate to these asymptotes would require a matrix size which is intractable and computationally very demanding to perform exact quantum dynamics simulations. Instead of including the full set of 10 electronic states correlating to the C(¹D) asymptote, we have only considered the 4 states corresponding to the two Δ components. This truncation will yield an incomplete description of some asymptotic states as we will discuss later.

Considering that the Λ quantum number and charge of each diatomic fragment are conserved across the diabatic states, a correlation between reactants and products can be done. On top of that, it is necessary to use three-body information to complete the correlation of states with the same Λ that tend to the same asymptote. The correlation diagram of the diatomic states is presented in Fig. 2.

The diabatic states are ordered as follows. The first two states are Σ and are different in the A′ and A″ diabatic manifolds. The next four states are of Π character and are degenerate in the A′ and A″ manifolds, correlating to the same products. The same occurs for the last two states, which are of Δ character. In the following we will use A′ and A″ to refer to adiabatic states and D′ and D″ for diabatic states. The states included in the model are enumerated in Table 1, together with the energies in the different minima and total dissociation.

Table 1 Minimum energies of the different diatomics and total dissociation channels included in the PEM. The energies are referred to the CH₂⁺ minimum. Each line represents a diabatic state in the PEM. When no number is indicated the electronic state does not have a minimum

State	[CH + H]^+	ΔE/eV	[C + H2]^+	ΔE/eV	[H + H + C]^+	ΔE/eV
1D′	CH^+(^1Σ^+) + H(^2S)	4.71	H2(^3Σu^+) + C^+(^2P)	—	H + H + C^+(^2P)	8.87
2D′	CH^+(^3Σ^+) + H(^2S)	—	H2(^1Σg^+) + C^+(^2P)	4.25	H + H + C^+(^2P)	8.87
1D″	CH(^2Σ^−) + H^+	11.17	H2^+(^2Σu^+) + C(^3P)	—	H^+ + H + C(^3P)	11.33
2D″	CH^+(^3Σ^−) + H(^2S)	9.49	H2^+(^2Σg^+) + C(^3P)	8.57	H^+ + H + C(^3P)	11.33
3D′, 3D″	CH^+(^3Π) + H(^2S)	5.90	H2(^3Σu^+) + C^+(^2P)	—	H + H + C^+(^2P)	8.87
4D′, 4D″	CH^+(^1Π) + H(^2S)	7.74	H2(^1Σg^+) + C^+(^2P)	4.25	H + H + C^+(^2P)	8.87
5D′, 5D″	CH(^2Π) + H^+	7.88	H2^+(^2Σu^+) + C(^3P)	—	H^+ + H + C(^3P)	11.33
6D′, 6D″	CH^+(^3Π) + H(^2S)	—	H2^+(^2Σg^+) + C(^3P)	8.57	H^+ + H + C(^3P)	11.33
7D′, 7D″	CH(^2Δ) + H^+	10.86	H2^+(^2Σu^+) + C(^1D)	—	H^+ + H + C(^1D)	12.68
8D′, 8D″	CH^+(^1Δ) + H(^2S)	11.33	H2^+(^2Σg^+) + C(^1D)	9.94	H^+ + H + C(^1D)	12.68

---

2.3 Potential energy matrix

Two independent PEMs are considered: for the A′ and A″ irreducible representations. The diagonal elements correspond to the diabatic states and the non-diagonal elements to their couplings. The Π and Δ subblocks are identical in both representations, as well as Π–Δ couplings, so the differences between the energies are only due to the Σ states and how they couple with the Π and Δ states in each representation.

In our model we describe 16 electronic states, 8 in each A′ and A″ representations. The PEM is factored into a zeroth order diagonal matrix (U⁰), which accurately includes all the possible diatomic and total dissociation states, plus a three-body neural network matrix to incorporate corrections to the diabatic states and their couplings in the 3-body region (U^NN),

U = U⁰ + U^NN (1)

The elements of U⁰ are expressed as:

U⁰_ii = V_1B + V_2B + V_lr (2)

where V_1B are the atomic energies of the total dissociation asymptote, V_2B is the sum of the three possible diatomics in the diabatic state and correspond to diabatized states of the CH⁺, CH, H₂⁺ or H₂ systems. V_lr is the long range interaction of the rearrangement channel.

The U^NN matrix elements are computed from a permutationally invariant polynomial neural network (PIP–NN).³⁸ Those associated to Π and Δ states are shared between the two PEM in order to maintain the degeneracy in linear configurations. The output of the neural network is multiplied by a damping function of the CH distances to make them zero as the system tends to a 2-body or total dissociation region. Additional constraints to the non-diagonal elements arise from the fact that for highly symmetric configurations (C_∞v, D_∞h) the PEM should resemble the adequate block structure. To impose these constraints, couplings between Σ–Π, Σ–Δ and Π–Δ are multiplied by a sine function of the H–C–H angle. The full diabatic matrices and further details on the neural network implementation are presented in the SI.

While in theory this method can be generalized to more than 3 particles, the main difficulty resides in the definition of U⁰. As the number of atoms increases, the number of relevant asymptotic states can be very large, which would lead to an intractable diabatic matrix. On the other hand, U^NN can be easily generalized to other systems with a PIP representation.^39–41

U^NN is trained by minimizing the following loss function ([script L]) which, in contrast to other methods, only requires adiabatic information and hence makes the diabatization method automatic:

where V_s and V_s,pred are the ab initio and predicted adiabatic energies for the electronic state s. The second term regularizes the training process forcing the non-diagonal terms to be small and avoiding too large couplings where not needed. In this case λ₂ = 10⁻³ so that non-diagonal couplings are only minimized if they do not affect the accuracy of the adiabatic energies. N is the total number of ab initio energies and n_a and n_d the number of adiabatic and diabatic states which are trained in each representation, respectively.

Two key aspects make this method different from similar PEM factorizations such as the PM-DDNN²⁶ and turn the diabatization method into fully automatic: (1) sharing the degenerate diabatic states and couplings among the A′ and A″ irreducible representations and (2) constraining the non-diagonal couplings at highly symmetric configurations to resemble the adequate Σ, Π and Δ subblock structure in linear configurations.

2.4 Features of the PEM

The PEM has been trained with n_a = 5 and n_d = 8 in eqn (3), meaning that the neural network term will modify the 8 diabatic states in each representation to finally reproduce accurately the 5 lowest adiabatic states of each representation. These states cover an energy up to 14 eV in the Frank–Condon region of the CH₂⁺, which is the relevant energy range for astrochemical considerations. In the diatomic and total dissociation regions the number of accurate states is in general larger than five thanks to the U⁰ term. The ab initio calculations have been performed with the state averaged CASSCF level of theory followed by a MRCI calculation with the cc-pCVTZ basis set with ORCA 6.0.^42–45 In all the calculations the core orbitals of the C atom are closed. A total of seven doublet states are computed for each representation of the C_s point group.

U⁰ is responsible of reproducing the full PEM up to a rather small distance between the three atoms as depicted in Fig. 3, where a comparison between U⁰ and the ab initio data is presented. For R = 6 Å we find a good agreement with the exception of the states that are not included in the zeroth order model: some that correlate with C(¹D) or the curve corresponding to H₂ + C⁺(⁴P) which is visible in the right panels. As the third atom approaches to distances ≈3 Å, U⁰ is still rather accurate and only for very short distances the zeroth order model is clearly in disagreement. It is in this short range region where U^NN acts, correcting the U⁰ diabatic energies and including the missing couplings.

Fig. 3 Black dots represent the ab initio energy and solid lines the diabatic energies predicted by U⁰. From top to bottom, represents the change of the diatomic curves as the third atoms approaches. In all cases θ = 140° and γ = 90°.

Fig. 4 displays a path connecting the CH₂⁺ in its linear configuration with the [CH + H]⁺ asymptote on the left and the [C + H₂]⁺ asymptote in the right. In the top panel the adiabatic energies are compared with the calculated ab initio, represented with dots. For the 5 lowest adiabatic states we find a very good agreement along the path. In the bottom panel, the evolution of the diabatic states towards the photofragments is depicted. For instance, the minimum in the CH₂⁺ configuration comes from the 5D′, 5D″ diabatic states which connect CH(²Π) + H⁺ and H₂⁺(²Σ_u⁺) + C(³P) asymptotes, in accordance with our previous correlation diagram for this system.²⁷ However, this is not evident from Fig. 3. In U⁰ the lowest diabatic Π state close to the CH₂⁺ minimum are the 3D′, 3D″ states, which correlate with CH⁺(³Π) + H(²S) and H₂(³Σ_u⁺) + C⁺(²P), so one could assume that U^NN will try to preserve this order. Instead, U^NN pushes the 3D′, 3D″ states up in energy while bringing down the 5D′, 5D″ to reproduce the minimum of the CH₂⁺. Further comparisons of the adiabatic and ab initio results can be found in the SI.

Fig. 4 Representation of a path joining the [CH + H]⁺ and [C + H₂]⁺ asymptotes from the CH₂⁺ linear configuration in the center of the plot. Towards the [CH + H]⁺ asymptote a C–H distances is elongated. Towards the [C + H₂]⁺ asymptote the extraction of the C between both H atoms is coordinated with a progressive decrease of the H–H distance towards its equilibrium geometry in the H₂ molecule. In the top and bottom panels the adiabatic and diabatic energies are presented, respectively. The dots in the top panel represent the ab initio energies for those geometries.

The root mean squared errors (RMSE) of the different states are presented in Table 2. The largest error in the fifth state is due to the missing Σ and Π components of the C(¹D) + H(²S) + H⁺ asymptote, which limits the capability of the model to describe the H₂⁺(²Σ_g⁺) + C(¹D) asymptote. Excluding [C + H₂]⁺ geometries, the errors of the fifth adiabatic state up to 14 eV are 86.0 meV and 76.4 meV for the A′ and A″ irreducible representations, respectively.

Table 2 RMSE for the five lowest adiabatic states in each representation. The errors are given in meV and in parenthesis the number of points with energy up to 14 eV with respect to the CH₂⁺ minimum, which have been used to compute the error

Irrep.	1	2	3	4	5
A′	35.0 (20 542)	41.3 (20 209)	46.1 (19 557)	76.4 (18 763)	523.8 (18 331)
A″	45.4 (20 421)	43.7 (19 963)	41.2 (19 239)	57.0 (18 500)	519.2 (17 682)

---

More relevant to the quantum dynamics are the non-adiabatic coupling matrix elements (NACMEs) which can be analytically computed from the PEM and compared with ab initio calculations in regions close to a CI.⁴⁶ Fig. 5 presents a comparison between those computed from the PEM and ab initio at a CASSCF level close to the CI between the first two adiabatic states in the CH⁺ + H exit.

Fig. 5 Non-adiabatic coupling matrix elements computed close to the CI (θ = 0.01°) between the first two adiabatic states in the CH⁺ + H exit. In blue those from the PEM and in red the ab initio computed at a CASSCF level. ξ represents the coordinate with respect to which the NACME is computed.

3 Photodissociation dynamics

The photodissociation dynamics is studied with the quantum time-dependent MADWAVE3 code⁴⁷ as described previously^48–51 using the CH + H Jacobi coordinates, with r being the HC internuclear vector, R joining HC center-of-mass to the second H atom, and γ = r·R/rR. The rovibrational H₂ products are analyzed using the reactant-product transformation method.⁵² The dynamics is performed in the 8 × 8 diabatic representations of the A′ and A″ symmetries independently.

Here we consider the absorption from the ground vibrational states in the adiabatic X̃²A′ and òA″ states, described in the SI, corresponding to initial total angular momentum J_i = 0. The initial wave packet is built by multiplying the bound state by the transition electric dipole moment as described previously.^48–51 The transition dipole moments have been calculated from the ground adiabatic states to the 5 (4) adiabatic excited electronic A′ (A″), as described in the SI. These adiabatic transition dipole moments are then transformed to the diabatic representation using the eigen-vectors obtained after diagonalization of the 8 × 8 diabatic matrix at each nuclear configuration. Thus, the simultaneous excitation of the first 5 (4) A′ (A″) adiabatic states is considered here.

We study the J_i = 0 → J = 1⁻ rotational excitation for 4 different electronic transitions: from either X̃²A′(0, 0, 0) or òA″(0, 0, 0) initial states towards 8 electronically excited and coupled diabatic states, of either A′ or A″ symmetry.

The absorption spectra obtained from the autocorrelation function^50,51 shows an excellent agreement with the total photodissociation flux over all fragments above the first dissociation limit, as displayed in Fig. 6. This serves not only to check to convergence of the calculation, but also to guarantee that all the flux over different fragmentation channels is properly collected.

Fig. 6 Photodissociation cross sections for the J_i = 0 → J = 1⁻ transition from the X̃²A′(0, 0, 0) (solid lines) panels) or òA″(0, 0, 0) (right panels) to the 8 × 8 A′ (bottom panels) or A″ (dashed lines) coupled diabatic states.

The assignment of the bands is done by comparing these spectra to those obtained in the adiabatic representation, as described in the SI. The most intense peaks are those arriving to 5A′ and 4A″ electronic states, because these excited states have considerably larger transition dipole moments with X̃²A′ and òA″ states than the other electronic states at lower energies.

The spectra to the A′ manifold of states obtained from the X̃ or the à states show significant differences, attributed to the different initial wave packets, formed from different initial bound states (specially in the angular coordinate) and different transition dipole moments (see SI).

The transitions towards the A″ states show a drop around 10 eV, which is clearly attributed to the energy separation among the A″ states, less densely packed than the A′ states. The drop at 10 eV observed in the òA″ → A′ transition is due to the reduction of the 1A″ → 3 and 4 A′ matrix elements of the dipole moments.

There is a great variety of fragmentation processes, as illustrated by the electronic states of the fragments shown in Fig. 2: total fragmentation

CH₂⁺(X̃²A′, òA″) + hν → C⁺ + H + H (4)

→ C + H⁺ + H,

formation of CH or CH⁺

CH₂⁺(X̃²A′, òA″) + hν → CH⁺ + H (5)

→ CH + H⁺,

and formation of H₂ or H₂⁺

CH₂⁺(X̃²A′, òA″) + hν → H₂ + C⁺ (6)

→ H₂⁺ + C.

These different fragmentation routes are distinguished by the analysis of different fluxes on each electronic state and on individual rovibrational states of H₂/H₂⁺ and CH/CH⁺ diatomic products.⁵² The total flux, obtained for the CH internuclear distance r = 6 Å, collects the total fragmentation, eqn (4), the formation of H₂ or H₂⁺, eqn (6), depending on the electronic state, and one half of CH or CH⁺, eqn (5), expressed as

The fluxes for H₂/H₂⁺ and one of the two CH/CH⁺ rearrangement channels are obtained by evaluating the sum over all possible rovibrational diatomic eigen states in each channel for each electronic states, giving the quantities F_H2+C and F_CH+H/2. For the [CH + H]⁺ there are two equivalent rearrangement channels giving the same flux, F_CH+H/2 each. In Fig. 7 the cross sections for each channel and final electronic states obtained after photodissociation of the X̃²A′(0, 0, 0) ground state towards the 8 lower A′ and A″ final electronic states, respectively, are shown with a similar qualitative behavior to those obtained from òA″(0, 0, 0) initial state shown in the SI.

Fig. 7 Photodissociation cross sections for the X̃²A′(0, 0, 0) → all D′ states and all D″ states, as indicated in each panel, separating individual fluxes for the different fragments of each rearrangement channel in each electronic state, as discussed in the text.

The main fluxes in the higher energy interval, E > 11 eV, are either to total fragmentation or to CH/CH⁺ diatomic fragments, with small proportion of H₂/H₂⁺ products. For the 5, 6, 7 and 8 diabatic D′ or D″ states, the total fragmentation yield neutral carbon atoms, in either ³P or ¹D states and the charge is in the hydrogen atom. On the contrary, for the 3, 4 D′ or D″ states the charge is in the carbon atoms in the C⁺(²P) state. The Σ states in the A′ and A″ symmetry show different behavior: in the A′ manifold the charge is in the C⁺(²P) ion, while in the A″ manifold the charge is in the hydrogen atoms, leading neutral C(³P) atoms.

The most astonishing result is the high cross section to form neutral CH products, in either the ground X²Π (states 5D′ and 5D″) or A²Δ (states 7D′ and 7D″). CH can be formed in the neutral–neutral reaction C(³P) + H₂, but is endothermic by ≈1 eV. Thus the neutral reaction can only occur when H₂ is initially vibrationally excited in v = 2 or 3. Abundant highly vibrationally excited H₂ is only found in strongly UV illuminated regions such as the borders of molecular clouds (photodissociation or PDR regions) or protoplanetary discs (PPDs). In such environments, the high UV flux will also induce the formation of neutral CH fragments by the photodissociation of CH₂⁺.

The formation of CH_n⁺ is attributed to the hydrogenation chain of reactions CH_n⁺ + H₂ → CH_n+1⁺ + H, which stops in CH₃⁺, which in turn reacts with many other molecules, being the origin of several chemical networks for the formation of many larger hydrocarbons. CH⁺ is one of the first observed molecules in the ISM³² and since then widely observed. CH₃⁺ was recently observed³³ in a protoplanetary disk (d203–506) illuminated by the strong far ultraviolet (FUV) radiation field from nearby massive stars in Orion's trapezium cluster. However, CH₂⁺ has not been observed so far.

Thus the confirmation of this hydrogenation chain requires the detection of CH₂⁺, and its abundance determination requires incorporating all the destruction pathways, such as photodissociation, dissociative recombination with electrons and hydrogenation.

The total photodissociation cross section of CH₂⁺(X̃) (summing the contributions towards A′ and A″ manifolds) is shown in Fig. 8 and compared to those of CH₃⁺ from ref. 36 and CH⁺, which is recalculated in this work as described in the SI and is in agreement with that reported by Kirby et al.³⁵ The photodissociation cross section of CH₂⁺ is larger than that of CH₃⁺ but smaller than that of CH⁺.

Fig. 8 Total photodissociation cross section for the ground vibrational state of CH₂⁺ (X̃²A′) (summed over all A′ and A″ states), compared to those of CH₃⁺ from ref. 36, and CH⁺, recalculated here (see text) and very similar to that of Kirby et al.³⁵

The photodestruction is determined by the photodissociation rate obtained as the integral of the photodissociation cross section with the energy dependent radiation field. Using Draine's mean interstellar radiation field,⁵³ the photodissociation rate of CH₂⁺ is 6.45 × 10⁻¹¹ s⁻¹, while that of CH₃⁺ is³⁶ 6.83 × 10⁻¹² s⁻¹ and for CH⁺ is 2.89 × 10⁻¹⁰ s⁻¹. The rates obtained here for CH₂⁺ and CH⁺ are slightly lower than the values reported for the ISRF field in ref. 3 of 1.4 and 3.3 × 10⁻¹⁰ s⁻¹, respectively. One reason is that the ISRF flux is a modification of the Draine's field. The difference for CH₂⁺ is more important because a simple vertical transition was assumed to calculate the photodissociation cross section.

The total photodissociation cross section to form neutral CH products is shown in Fig. 9, which are rather high specially for CH(X²Π) at higher photon energies, where the UV radiation field is not shielded. The photodissociation rates to form CH in the ²Π and ²Δ electronic states are 9.96 × 10⁻¹² s⁻¹ and 3.37 × 10⁻¹² s⁻¹, respectively, using Draine's mean interstellar radiation field.⁵³ These values are small, but considering that the H₂ + C neutral reaction is endothermic by 1 eV, it can be an alternative source of CH.

Fig. 9 Total photodissociation cross section to form neutral CH products from the ground vibrational state of CH₂⁺ (X̃²A′ (summed over all A′ and A″ states).

4 Conclusions

The study of the branching ratios among the different photoproducts in photodissociation over large energy intervals, as needed in astrochemistry, requires the description of non-adiabatic dynamics, involving many electronic states. To overcome the singularities of non-adiabatic couplings at conical intersections, the most widely used method is a unitary transformation to a diabatic representation, which is not unique. There are many diabatization methods, but usually they are restricted to few states. In this work we have developed a diabatization method describing up to 16 electronic states, divided in two symmetry blocks. The method has been applied to the photodissociation of CH₂⁺, not yet dynamically studied, and that present many photoproducts, bound and dissociative and with several charge distributions.

The diabatization method consists in building a basis of electronic states to represent the possible products, with different charges and electronic excitation. A correlation among the rearrangement channels is done, imposing symmetry restrictions, such as the conservation of the projection of angular momentum, Λ, on the diatomics fragments. The full diabatic matrices are then factorized in two terms, a diagonal zero-order term and a full matrix built using a neural network method, V_NN. The zero-order term allows to describe the diagonal interaction among the fragments up to intermediate distances, ≈3–4 Å, keeping the nature of the asymptotic states. The V_NN matrix is build with an artificial neural network, imposing symmetry restrictions on the couplings at symmetric configurations and introducing a penalty function to force non-diagonal terms to be as small as possible. This process allowed an accurate description of the lowest 10 adiabatic eigenvalues in the short distance range and 16 in the asymptotes, separated in two symmetry blocks, covering energies up to 14 eV.

The diabatic model allows the quantum dynamics of the total photodissociation cross sections rate constants for CH₂⁺ up to 13.6 eV, interval needed to astrochemical models. More important, it allows to describe the partial cross section for a great variety of fragmentation processes, with the products in different arrangement channels, charge and electronic states: total fragmentation, C⁺ + H + H and C + H + H⁺, CH⁺(X¹Σ,) + H, CH(X²Π, ²Δ), H₂(X¹Σ,) + C⁺ and H₂⁺(²Σ_g).

Here we focus in the formation of the CH radical, whose formation in the C + H₂ reaction is endothermic. In strongly illuminated astronomical objects, the CH₂⁺ photodissociation leads to significant of CH radical, and it is important to determine its contribution in PDRs and PPDs under different conditions to determine its contribution.

CH₂⁺ is a prototype and this kind of studies are needed to search for alternative formation routes of many molecules and radicals in astronomical objects illuminated by UV radiation. For this purpose diabatization is a key important requirement, specially in cases where the fragments present crossings. At these crossings at long distance the non-adiabatic couplings vanish and surface hopping methods would adiabatically continue in one adiabatic state, corresponding to a mixture of two states each one assigned to one chemical product, as illustrated in Fig. 2 for the present case.

Author contributions

Pablo del Mazo-Sevillano: conceptualization (equal); methodology (equal); software (equal); writing – review & editing (equal). Alfredo Aguado: conceptualization (equal); funding acquisition (equal); methodology (equal); writing – review & editing (equal). Susana Gómez-Carrasco: conceptualization (equal); funding acquisition (equal); writing – review & editing (equal). Octavio Roncero: conceptualization (equal); funding acquisition (equal); methodology (equal); project administration (equal); software (equal); writing – review & editing (equal).

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp01221c.

The total photodissociation cross sections calculated here for CH₂⁺ and CH⁺ are available at ZENODO https://zenodo.org/records/19847742https://zenodo.org/records/19847742. The codes used for the calculations have been properly cited, and the details of the calculations are provided in the supplementary information (SI). Any other data that support the findings of this study are available from the authors, upon reasonable request.

Acknowledgements

This work has received funding from Ministerio de Ciencia, Innovación y Universidades, MICIU (Spain), under Grants No. PID2021-122549NB-C21, PID2021-122549NB-C22 and PID2024-156686NB-I00. Computational assistance was provided by the Supercomputer facilities of Lusitania founded by the CénitS and Computaex Foundation.

Notes and references

1. M. C. van Hemert and E. F. van Dishoeck, Chem. Phys., 2008, 343, 292.
2. E. F. van Dishoeck and R. Visser, Laboratory Astrophysics: from Molecules through Nanoparticles to Grains, Wiley-VCH, 2015, p. 1.
3. A. N. Heays, A. D. Bosman and E. F. van Dishoeck, Astron. Astrophys., 2017, 602, A105.
4. H. R. Hrodmarsson and E. F. van Dishoeck, Astron. Astrophys., 2023, 675, A25.
5. D. R. Yarkony, Modern electronic structure theory, World Scientific, Singapore, 1995, p. 642.
6. D. R. Yarkony, Conical Intersections: electronic structure, dynamics and spectroscopy, Advanced series in Physical Chemistry, World Scientific Publishing Co., 2004, p. 175.
7. F. T. Smith, Phys. Rev., 1969, 179, 111.
8. K. Ruedenberg and G. J. Atchity, J. Chem. Phys., 1993, 99, 3799.
9. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1982, 77, 6090.
10. H. Köppel, Conical Intersections: electronic structure, dynamics and spectroscopy, Advanced series in Physical Chemistry, World Scientific Publishing Co., 2004, p. 175.
11. A. Thiel and H. Köppel, J. Chem. Phys., 1999, 110, 9371.
12. T. Pacher, H. Köppel and L. S. Cederbaum, J. Chem. Phys., 1991, 95, 6668.
13. T. Pacher, L. S. Cederbaum and H. Köppel, Adv. Chem. Phys., 1993, 84, 293.
14. B. N. Papas, M. S. Schuurman and D. R. Yarkony, J. Chem. Phys., 2008, 129, 124104.
15. S. Ghosh, S. Mukherjee, B. Mukherjee, S. Mandal, R. Sharma, P. Chaudhury and S. Adhikari, J. Chem. Phys., 2017, 147, 074105.
16. S. Mukherjee, S. Ravi, K. Naskar, S. Sardar and S. Adhikari, J. Chem. Phys., 2021, 154, 094306.
17. S. Hazra, S. Mukherjee, S. Ravi, S. Sardar and S. Adhikari, ChemPhysChem, 2022, 23, e202200482.
18. Y. Guan, H. Guo and D. R. Yarkony, J. Chem. Phys., 2019, 150, 214101.
19. Y. Guan, D. H. Zhang, H. Guo and D. R. Yarkony, Phys. Chem. Chem. Phys., 2019, 21, 14205–14213.
20. Y. Guan, C. Xie, D. R. Yarkony and H. Guo, Phys. Chem. Chem. Phys., 2021, 23, 24962–24983.
21. D. M. G. Williams and W. Eisfeld, J. Chem. Phys., 2018, 149, 204106.
22. C. Xie, X. Zhu, D. R. Yarkony and H. Guo, J. Chem. Phys., 2018, 149, 144107.
23. W. Li, Y. Liang, X. Niu, D. He, W. Xing and Y. Zhang, J. Chem. Phys., 2024, 161, 044310.
24. W. Li, B. Dong, X. Niu, M. Wang and Y. Zhang, J. Chem. Phys., 2024, 074302.
25. Y. Shu and D. G. Truhlar, J. Chem. Theory Comput., 2020, 16, 6456.
26. Y. Shu, F. B. Akher, H. Guo and D. G. Truhlar, J. Phys. Chem. A, 2024, 128, 1207–1217.
27. P. del Mazo-Sevillano, A. Aguado, F. Lique, R. A. Jara-Toro and O. Roncero, Phys. Chem. Chem. Phys., 2025, 27, 15775–15786.
28. S. Sršeň, O. A. von Lilienfeld and P. Slavíček, Phys. Chem. Chem. Phys., 2024, 26, 4306.
29. P. Jensen, M. Brumm, W. P. Kraemer and P. R. Bunker, J. Mol. Spectrosc., 1995, 172, 194.
30. J. H. Black and A. Dalgarno, Astrophys. J., 1976, 203, 132–142.
31. D. Smith, Chem. Rev., 1992, 92, 1473.
32. A. E. Douglas and G. Herzberg, Astrophys. J., 1941, 94, 381.
33. O. Berné, M.-A. Martin-Drumel and I. Schroetter, et al., Nature, 2023, 621, 56.
34. R. P. Saxon, K. Kirby and B. Liu, J. Chem. Phys., 1980, 73, 1873.
35. K. Kirby, W. G. Roberge, R. P. Saxon and B. Liu, Astrophys. J., 1980, 855.
36. P. del Mazo-Sevillano, A. Aguado, J. R. Goicoechea and O. Roncero, J. Chem. Phys., 2024, 160, 184307.
37. G. Theodorakopoulos and I. D. Petsalakis, J. Mol. Struct.:THEOCHEM, 1991, 230, 205.
38. B. Jiang and H. Guo, J. Chem. Phys., 2013, 139, 054112.
39. J. Li, Z. Varga, D. G. Truhlar and H. Guo, J. Chem. Theory Comput., 2020, 16, 4822–4832.
40. P. del Mazo-Sevillano, D. Félix-González, A. Aguado, C. Sanz-Sanz, D. Kwon and O. Roncero, Mol. Phys., 2024, e2183071.
41. P. L. Houston, C. Qu, Q. Yu, P. Pandey, R. Conte, A. Nandi and J. M. Bowman, J. Chem. Theory Comput., 2024, 20, 3008–3018.
42. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2025, 15, e70019.
43. C. Kollmar, K. Sivalingam, B. Helmich-Paris, C. Angeli and F. Neese, J. Comput. Chem., 2019, 40, 1463–1470.
44. F. Neese, J. Comput. Chem., 2022, 44, 381.
45. M. Ugandi and M. Roemelt, Int. J. Quantum Chem., 2023, 123, e27045.
46. C. Sanz-Sanz, A. Aguado, O. Roncero and F. Naumkin, J. Chem. Phys., 2015, 143, 234303.
47. O. Roncero and P. del Mazo-Sevillano, Comput. Phys. Commun., 2025, 109471.
48. M. Paniagua, A. Aguado, M. Lara and O. Roncero, J. Chem. Phys., 1999, 111, 6712.
49. A. Aguado, M. Paniagua, C. Sanz-Sanz and O. Roncero, J. Chem. Phys., 2003, 119, 10088.
50. A. Chenel, O. Roncero, A. Aguado, M. Agúndez and J. Cernicharo, J. Chem. Phys., 2016, 144, 144306.
51. A. Aguado, O. Roncero, A. Zanchet, M. Agúndez and J. Cernicharo, Astrophys. J., 2017, 838, 33.
52. S. Gómez-Carrasco and O. Roncero, J. Chem. Phys., 2006, 125, 054102.
53. B. T. Draine, Astrophys. J., Suppl. Ser., 1978, 36, 595.

---