A time-independent, variational method for studying the photodissociation of triatomic molecules

Abstract

The photodissociation of molecules is becoming an increasingly important factor to consider in the evolution of exoplanets' atmospheres orbiting around UV-rich stars, as it leads to the enrichment of atmospheric complexity. A new method is developed for computing the rotationally and vibrationally resolved photodissociation spectrum of triatomic molecules. The time-independent Schrödinger equation is solved using the variational nuclear motion program EVEREST; a new code EXOCSMOOTH is employed to compute the cross-sections by applying Gaussian smoothing to a set of discrete transitions into the continuum. HCN is chosen as the test molecule, as it has been widely studied in the literature. Results are compared with the available experiments. Temperature dependence is explored for temperatures up to 2000 K.

---

1 Introduction

Understanding the photochemistry of planetary atmospheres is becoming extremely important in this new revolutionary era of astronomy. The recent discovery of photochemically produced SO₂ in WASP-39b is an example where photochemical processes, including photodissociation, play a crucial role in explaining the complexity of atmosphere composition and dynamics, as well as non-LTE effects.^1,2 Exoplanets orbiting near M-dwarfs, small oxygen-rich stars with temperatures between 2300 and 3800 K which are the commonest stars in our galaxy, experience high UV fluxes of photons from their host stars.^3–6 Due to the planets’ proximity to the star, they have temperatures on the order of magnitude of thousands Kelvin, generating vibrationally and rotationally excited molecules. Challenges arise in collecting photodissociation cross-sections at these high temperatures,⁷ due to the lack of experimental and computational data. Experiments struggle to reach the temperatures required to model the top of hot atmospheres, and calculations often assume a harmonic ground state wave function.⁸ As a reflection of these two factors, standard databases contain data for molecules at interstellar temperatures (close to 0 K).^9,10

In recent years, cross-section calculations have begun to include rovibrational excitation of molecules as a function of temperature.¹¹ The ExoMol database^12,13 was designed to produce comprehensive line lists of hot bound–bound transitions for molecules present in exoplanets’ atmospheres and is expanding to include photodissociation data. Previous studies¹⁴ have demonstrated that the photodissociation of diatomic molecules can be simulated using the program Duo,¹⁵ which was originally developed to solve the bound–bound nuclear motion problem for diatomics allowing for different types of couplings and crossings between potential energy curves. Duo uses a grid-based variational solution of the Schrödinger equation based on the use of a sinc-DVR basis which has been shown to be highly accurate¹⁶ and which results in a discretised representation of the continuum. ExoCross¹⁷ is used to generate the raw cross-section for different grid sizes. The results obtained using different radial grids and smoothing the results with an appropriate Gaussian function generate spectra comparable with those of experiments. This methodology was successfully applied to produce photodissociation line lists for HCl and HF;¹⁸ these results were subsequently corroborated by Qin et al.¹⁹

The transition from two to three atoms adds complexity to the computations involved in generating photodissociation cross-sections, not least because the continuum now becomes multi-dimensional. Determining the best set of internal coordinates for the given potential energy surfaces and defining the transition dipole moments are non-trivial tasks that depend on the chemical species under examination. HCN was selected as the initial molecule to assess the generation of photodissociation cross-sections for triatomic molecules as at lower energies it has a clear set of photodissociation products: H + CN. HCN has an ubiquitous presence in space; it is detected in the interstellar medium,^20–24 comets,^25–27 planets,^28–30 exoplanets³¹ and, finally, in the atmosphere of Titan,³² where it is supposed to be a main reservoir of carbon.^33,34 The HCN/HNC ratio is very sensitive to stellar UV flux^35,36 and can be used as a tracer for the heating of the nebular molecular gas by UV photons.^37–39

The photochemistry of HCN has undergone extensive examination, encompassing both experimental and computational approaches. Experiments explored dissociation from the predissociative 2 ¹A′ ← 1 ¹A′ transitions at various wavelengths, successfully reconstructing the potential dissociation mechanism.^40–43 Dissociation events from 1 ¹A′′, 2 ¹A′′, and 2 ¹A′ states involve the cleavage of the H–CN bond through bending progression.⁴⁴ Analysis of the electronic state distribution of the CN fragment, originating from the bending of the 3 A′ ← 1 A′ transition, indicates potential intersystem crossing and/or internal conversion between various excited states.⁴⁵ This hypothesis has been corroborated by measurements from ¹A′(v_CH = 3) to the 1 ¹A′′ electronic states: the lack of rotational population suggests that the bond breaks when HCN is in the linear arrangement.⁴⁶

The ground state potential energy surface (PES) has been computed by various research groups,^47–52 whereas surfaces from the excited states are relatively scarce. The first study of excited states was published by Perić et al.⁵³ in 1987. This work was followed by studies performed at the University of New Mexico between 2001 and 2003 on predissociation from the 1 ¹A′′ (ref. 54–56) and 2 ¹A′ (ref. 57) electronic states. Vibrational predissociation resonances in the latter state accelerate the dissociation due to efficient energy flow facilitated by an accidental coincidence between H–C and C–N vibrational normal modes. Nayak and collaborators calculated a variety of excited singlet states,⁵⁸ followed by the work from Chenel et al.,⁵⁹ which was used to explore the rates in the interstellar medium.⁶⁰ The triplets states have been recently observed by Priyadarshini and coworkers.⁶¹

An alternative method to generate photodissociation cross-sections is through the use of the time-dependent Schrödinger equation. In this approach, the cross-section is proportional to the Fourier transform of the time-dependent autocorrelation function of the time evolution of the wave packet in an excited electronic state. Among all time-dependent methodologies, the multi-configurational time-dependent Hartree (MCTDH) approach, developed by Meyer et al.,⁶² stands out as one of the most versatile methods, as demonstrated by multiple improvements to the original code and its various implementations.^63,64 The first photodissociation application appeared, for NOCl⁶⁵ and NO₂,⁶⁶ two years after its introduction. These studies were quickly followed by investigations of the photodissociation of larger molecules, such as methyl iodide, which involved five nuclear degrees of freedom across three excited electronic states.^67,68 The recent perspective from Han et al.⁶⁹ discusses the most up-to-date applications of the method.

Although acknowledging the advantages of time-dependent methods for studying photodissociation, we opt for a time-independent approach for three reasons: (1) the ability to generate state-to-state photodissociation spectra, which are useful for modelling non-LTE effects, (2) the convenience of applying it to indirect photodissociation processes, as it is less time-intensive than the time-dependent approach⁷⁰ and (3) it allows the use of a single, often complicated, model for all photon absorption processes,⁷¹ which includes the calculation of photoabsorption cross-sections across an extended wavelength range.^72,73

In this work, we explore the ability to extend the approach used for diatomics to study continuum processes in larger systems. The DVR3D⁷⁴ and EVEREST⁷⁵ codes were chosen because they have been successfully used to describe bound–bound electronic transitions in previous studies, as demonstrated by Zak et al.⁷⁶ and Owens et al.,⁷⁷ respectively.

This paper is organised as follows: in the Methods section, a description and tests of our approach, the potential energy surface used, the choice of the basis, and the cross-section calculation are discussed. In the Results section, we compared our cross-section and the integrated intensities with previous experiments and calculations. Finally, the temperature dependence of photoabsorption is explored for three different temperatures.

2 Methods

2.1 Potential energy and transition dipole moment surfaces

Potential energy surfaces for this work are taken from the literature. The ground state has been taken from the study by Makhnev et al.⁵² and chosen for its spectroscopic accuracy when compared to experiments. The PESs and transition dipole moment surfaces (TDMS) for the excited A′ states have been provided through private communication with Chenel et al. (16ChRoAg).⁵⁹ The potential energy surfaces are represented using the bond-length-bond-angle (BLBA) representation with the following the convention: r₁ = r_CH, r₂ = r_CN and θ = HĈN. The transition dipole moments are oriented such that the z-axis is aligned to r_CN, the x-axis is perpendicular to the z-axis inside the molecular plane, and the y-axis component, which has a zero contribution, is perpendicular to the molecular plane. Fig. 1 shows the molecular coordinates and the reference axis of the transition dipole moment.

Fig. 1 BLBA coordinates of HCN are used in this work. The reference axis, shown below the molecule, indicates the transition dipole moment orientation.

Transitions from the 1 ¹A′ ground electronic state to the 2 ¹A′, 3 ¹A′, 4 ¹A′, and 5 ¹A′ excited states are considered. 2 ¹A′ and 3 ¹A′ states are mainly dissociative states with a bound region near the linear structure due to an energy barrier along the r_CH coordinate for both surfaces. No transitions between excited states are included because of their large excitation energies. The PESs are illustrated in Fig. 2. Details of these PESs are shown in Table 1.

Fig. 2 Cuts through the potential energy surfaces of the low-lying singlet electronic states of HCN considered in this work as a function of r_CH (up) and Θ (bottom).

Table 1 Linear configuration mapping, excitation energies (T₀) and vertical excitation energies (ΔE) for the HCN electronic states considered in this work. T₀ values in parenthesis are from 16ChRoAG⁵⁹

State	Linear configuration	T 0 (cm^−1)	ΔE (cm^−1)
1 ^1A′	X^1Σ^+	3418 (3468)	0
2 ^1A′	1^1Δ	54 040 (55 491)	70 241
3 ^1A′	1^1Π	65 916 (65 815)	75 162
4 ^1A′	2^1Π	75 744	82 534
5 ^1A′	3^1Π	85 389	86 203

---

2.2 Choice of the basis

Calculations of the photoabsorption cross-section were performed using DVR3D⁷⁴ and EVEREST.⁷⁵ The original DVR3D program suite, based on the use of the discrete variable representation (DVR) basis, was modified by Zak et al.⁷⁶ to calculate rovibronic transitions and tested for SO₂. The program suite EVEREST was used to solve the Schrödinger equation using a sinc-DVR basis.^78,79 The EVEREST input parameters were tuned to reproduce the energy levels published by Makhnev et al.⁵² for the J = 0 vibrational energies below 10 000 cm⁻¹, the two models agree with an RMSE (root-mean-squares-error) of 0.012 cm⁻¹.

For DVR3D, the basis set parameters used for the ground state and excited states must be the same; this is a necessary requirement for calculating the overlap between the two wave functions. DVR3D offers two forms of radial basis functions which can be optimised for a given system: Morse-like oscillators⁸⁰ and spherical oscillators;⁸¹ standard ground state HCN calculations such as those of Makhnev et al.⁵² use Morse-like oscillators for both radial coordinates. Here, to allow for dissociation, we experimented using spherical oscillators for the H–CN coordinate in Jacobi coordinates because they have better space filling properties. However, unlike the sinc-DVR basis which fill a finite, well defined region, the spherical oscillators in principle fill all space but the representation becomes increasingly sparse at large separations. Numerical experiments showed that, while it was possible to recover some of the photodissociation behaviours using spherical oscillators, this was less satisfactory than the results using sinc-DVR basis functions, which have the additional advantage that one can perform repeat calculations with shifted box sizes and obtain useful information.^71,82 The sinc-DVR basis successfully reproduced the energy levels and transitions observed in HCN, see below. The reason behind the different behaviour of the two basis functions is the way the tridimensional grid is built: the DVR basis set is constructed using a non-uniform, in principle semi-infinite Gaussian quadrature, while the sinc-DVR basis uses a finite and uniform grid. Fig. 3 shows how the distribution of the vibrational levels changes between a spherical oscillator DVR and a sinc-DVR basis for the 3 ¹A′ dissociative state. The two basis gave different results for energies above 70 000 cm⁻¹. The divergence is due to the spherical oscillator DVR basis being unconstrained in the radial coordinates: the basis explores larger regions which means that attempting convergence by increasing the number of points simply results in the production of an increasing, potentially infinite, number of states. For this reason, all calculations in this paper were performed using EVEREST.

Fig. 3 The first 600 vibrational energy levels of HCN (3 ¹A′). For energies below 70 000 cm⁻¹, spherical oscillator DVR and sinc-DVR basis results show a similar trend, but they diverge for higher energies.

2.3 Calculation of the photoabsorption cross-section

In this section, we provide a brief overview of the procedure and steps for generating cross-sections using EVEREST. A more detailed explanation of the code and underlying algorithms can be found in the original publication by Mitrushchenkov.⁷⁵

The first step required to calculate cross-sections is to generate the vibrational basis set for the electronic states of interest. This is done using the evvib.e routine. For each individual state, the corresponding symmetry is specified in terms of the electronic spin (S), the molecular term symbol (A), and the projection of the orbital angular momentum along the internuclear axis (Λ). The atomic compositions, molecular masses and coordinate systems are then specified. For each coordinate, the basis representation, the grid interval and the number of points need to be specified. The bond length coordinates are described using a sinc-DVR basis using the desired number of points (N^r_points), ranging from r_min to r_max. The bond angle is described by points using Legendre polynomials ranging from 0 to 180 degrees. This produces the vibrational levels for each k (projection of the angular momentum on the z-axis) of interest for a Hamiltonian of the desired size.⁸³

The second step consists of generating the rotational wave function using evrot.e. In this routine, the wave function for the desired J_s is generated. The rotational embedding is specified in this step.

In the third step, the dipole moment transition matrix is generated with the evdip.e routine, using the x and z-components of the TDMS. The routine evrli.e combines the wave function with the transition dipole moments. The evrs.e routine generates the molecular spectrum. The output is then transformed and stored in the ExoMol format.^13,84

A significant advantage of EVEREST is that potentials, dipole moments, and possible couplings are loaded as dynamic libraries, without the need to recompile the code for each calculation.

The EVEREST parameters used in the present study are reported in Table 2. The first section of the table presents the parameters for the vibrational input. Each state considered in this work is a singlet with A′ symmetry, while Λ varies depending on the state considered. Atom names and (atomics) masses are the input with the central atom last. Each calculation was made with a Hamiltonian of size 6000, specifying the first 100 roots for the ground state and the first 1000 roots for the electronically excited states. In the rotational input, z-embedding is specified with the z-axis lying along the r_CN bond. For both r_CH and r_CN, a radial grid spanning from 0.6 Å to 3.5 Å is used to take into account wavefunctions for highly excited vibrational states. At this step, all energy levels below 10⁴ cm⁻¹ for the ground state and below 10⁵ cm⁻¹ for the excited states are selected. The input runs through the first N J and K roots. The overall workflow is shown in Fig. 4.

Table 2 Input parameters used in the EVEREST code for HCN

Vibrational input
N pots	2
State 1 symmetry	1 A′ Λ1
State 2 symmetry	1 A′ Λ2
Atoms	H N C
Masses [Da]	1.007825
	14.003074
	12.000000
Geometry	BLBA
r 1 Grid type	Sinc-DVR
r 1 N points	50
r 1 Interval [Å]	0.6
	3.5
r 2 Grid type	Sinc-DVR
r 2 N points	50
r 2 Interval [Å]	0.6
	3.5
Θ Grid type	Legendre
Θ N points	50
Θ Interval [degrees]	0
	180
K min	0
K max	N
Hamiltonian size	6000
Number of roots	100
	1000
Rotational input
J interval	0-N
Embedding	z-Embedding
Embedding axis	r x
Diagonalisation	Jaco
Hamiltonian size	6000
E ^maxstate [ cm^−1]	10 000
	100 000
K min	0
k max	N
Dipole input
Number of dipole files	2
Alx	1
Embending	z-Embedding
Embending axis	r x

---

Fig. 4 Graphical representation of the workflow used in this paper. For each step, the program used in EVEREST is written in parenthesis in an italic character.

All calculations are performed in the UCL Department of Physics and Astronomy theory cluster using an Intel compiler and an OpenMP parallelisation scheme on 16 nodes.

2.4 Cross-section smoothing

With a sinc-DVR basis, one performs a calculation over a finite region (box) defined by internuclear distances r_min, and r_max which effectively places infinite potential barriers at the boundaries. For diatomic molecules,¹⁴ it is computationally cheap to repeat the calculation using a sinc-DVR basis generated with an increasing box size. This can form the basis of a stabilization method^85–87 and is being used by the ExoMol group to study predissociation.^71,73,82 The choice of the box size sets boundary conditions and controls the density resulting in a discretized continuum spectrum. By adjusting the size of the box by varying r_max, one can sample different areas of the continuum and produce cross-sections with a dependence on the calculation of the box size. For diatomics, it is computationally possible to perform hundreds or thousands of nuclear motion calculations⁷¹ as individual calculations are relatively quick, requiring, typically, less than an hour of computing time per iteration. This is not the case with triatomic molecules or larger polyatomics, where calculations can take days, weeks, or months per iteration, depending on the size and complexity of the model. As such, performing stabilization-style calculations is currently computationally prohibitive.

A practical alternative to the stabilization method is to compute a discretized continuum spectrum from a single-box calculation and smooth the corresponding temperature-dependent cross-sections, σ(ν), by applying a uniform Gaussian line profile with an appropriate half-width at half-maximum (HWHM, α). This is, at the moment, the only viable method for producing these cross-sections for polyatomics using EVEREST. The choice of α can be constrained by the structure of the discretized spectrum and the shape of the continuum cross-section.

Following Pezzella et al.,¹⁴ the initial cross-sections of 1 ¹A′→ U electronic transitions are produced with α = 10 cm⁻¹. Here, U represents any given unbound electronic state. These are the cross-sections which undergo smoothing.

To quantify optimum smoothness, we require that the resulting smoothed cross-section must be continuous and that the minimum value of α is used to prevent the over-distribution of intensity to the wings.

To assess the continuity of the cross-sections, one notes that, at their intensity maxima, the cross-sections have isolated turning points. In this region, the intensity distribution is bell shaped and, in general, direct photodissociation to purely repulsive states appears to have this shape as found by Pezzella et al.^14,18

As such, the smoothing of the cross-section can be constrained by the number of expected maxima, N_M. The continuity condition is so defined to be that the resulting, smoothed cross-section should have a number of roots in its first and second differential equal to N_M and N_M + 1, respectively, hence, should be twice differentiable.

It is found that there are usually a range of α values that meet this condition, and so the minimum value is taken to satisfy the condition on minimum smoothing. Eqn (1) expresses these conditions formally:

where N(f) is the number of roots as a function f. Due to the discrete nature of this problem, one cannot implement this procedure as a standard minimization problem through treatments like gradient descent and instead it must be solved iteratively.

The numerical method behind this is to establish a range of α values over which to test. At each iteration, the first and second forward derivatives are computed as follows:

where s is the sampling window which is chosen to deal with the numerical noise of producing these differentials. In general, a value of , where is the size of the wavenumber grid, is found to be appropriate for both first and second derivatives. Roots are counted by scanning the available smoothed cross-section and finding where the sign of the derivatives change. In areas of low intensity in the initial cross-section, small changes in intensity result in large changes of the gradient. While roots in these regions are not false positives, it is impractical to count them. To counteract this, a threshold value, K, is imposed such that a root will only be counted if the smoothed cross-section is [small sigma, Greek, circumflex](ν;α) >K.

Each iteration checks in the number of counted roots which meet the threshold conditions, upon the continuity condition being met or all test values are used. An example showing how the increasing α parameter affects the derivatives and shape of a toy model's cross-section is shown in Fig. 5. As can be seen, both α = 1000 cm⁻¹ and α = 10 000 cm⁻¹ meet the continuity condition imposed by eqn (1); however, the shape of the cross-section in the α = 10 000 cm⁻¹ case does not faithfully represent the shape of the raw cross-section in the top left panel; therefore, one must include the requirement for minimal smoothing. In this dataset, the threshold value was estimated as K = max(σ)/10.

Fig. 5 Toy model: Dependence of the shape of a photodissociation cross-section and its first and second derivatives on α in the smoothing procedure. Left panels show the smoothed cross-section, middle panels are the first derivatives, and right panels are the second derivatives.

The benefits of this procedure are that it is simple to implement, requires minimal computational expense, and provides good results in terms of smoothing. A limitation of this method is that upon meeting the aforementioned criteria, smoothing does not tend to be sufficient at higher ν. As ν increases, the gaps between successive vibrational bands also increases and so a minimal smoothing parameter as defined above would still allow for discrete structures to show through at high ν. One could increase α to account for this; however, this results in a slight over-smoothing in the lower ν region and is more challenging to automate reliably. This is a minor limitation, however, as the temperature dependent effects of these cross-sections will manifest themselves primarily in the lower wavenumber regions.

A Python3 code, ExoSmooth, has been written to automate the procedure of optimizing the HWHM. The program applies Gaussian profiles to the unsmoothed cross-section, starting at a given HWHM and increasing the HWHM by some user-defined amounts at each iteration until the continuity conditions described above are met. This ensures that the minimum amount of smoothing possible is applied to the cross-section. Indeed, as the Gaussian profile is normalised, the integral of the smoothed and unsmoothed cross-sections is equivalent; hence, photodissociation rates calculated with either cross-section are also equivalent. The documented code is made available from https://github.com/exomol. We note that the test performed as part of our initial study on diatomics¹⁴ showed that photodissociation rates obtained with the smoothed and unsmoothed data were the same, and that comparisons with (rotation-free) full photodissociation calculations also gave essentially the same photodissociation rate at low temperatures.

3 Results

3.1 Comparison with experiments

Of the four excited states available, there are no experiments available for the 2 A′ ← 1 A′ transition; the 3 A′ ← 1 A′ transition has been explored by different authors,^45,88,89 thanks to its well defined structure and being isolated from other electronic states; the 4 A′ ← 1 A′ and 5 A′ ← 1 A′ transitions are localised in the same energy region, making difficult to separate the two contributions in the absorption spectrum.^45,89,90

The 3 A′ is the only excited state which shows a clear vibrational progression. These levels make it possible to evaluate the quality of the potential used and the overall validity of the model. The quality of our calculations is tested against the computational results of 16ChRoAg and the experiments of Nuth et al. (82NuGl).⁸⁹ The first two models are computed including only the J = 1 ← 0 transitions, which effectively models T = 0 K. The results are shown in Fig. 6. The peak positions agree, in terms of relative intensities between the bending progression, with the measurements of 82NuGl, while the calculations of 16ChRoAg tend to be blue shifted with respect to the other sources. The cross-section magnitudes agree with 16ChRoAg up to v₂ = 6, while, for v₂ = 7,8 we observe a major splitting between the K components compared to the other two methods, leading to a lower cross-section peaks than 16ChRoAg. The splitting becomes evident for v₂ ≥ 9 in all sources. A clear drop of the cross-section is observed for energies higher than the 3 A′ dissociation barrier at 138.84 nm (72 025 cm⁻¹). The absorption spectrum presents both bound–bound features, due to the bending progression, together with an increasing dissociative character raising from dissociative threshold at wavelengths shorter than 138.84 nm.

Fig. 6 Comparison between our photoabsorption spectrum for the 3 A′ ← 1 A′ transition at T = 0 K and J = 0,1, with the computed model from 16ChRoAg⁵⁹ and the experiments from 82NuGl.⁸⁹ The blue dashed line separates the bound–bound spectrum at longer wavelengths from the photodissociation region at shorter wavelengths.

Table 3 reports measured vibrational band centers for the 3 A′(0,v₂,0) ← 1 A′(0,0,0) bands from MacPherson & Simons (78MaSi),⁸⁸ Lee (80Lee),⁴⁵ 82NuGl and our values at T = 0 K. The experimental vibrational values vary by 10 to 70 cm⁻¹ between the different measurements. In our calculations, the ground state transition is shifted by 46 cm⁻¹ from 80Lee with a RMSE of 68 cm⁻¹. The RMSE is 64 cm⁻¹ compared to that of 78Masi, while a lower RMSE of 47 cm⁻¹ is found when comparing the measurements of 82NuGl.

Table 3 Position for the peaks of the HCN (3¹A′ (0, v₂, 0) ← X̃′Σ⁺ (0,0,0)). Experimental values are from 78MaSi,⁸⁸ 80Lee⁴⁵ and 82NuGl.⁸⁹ Transition frequencies and the RMSE of the observed calculated results are in cm⁻¹

v 2	ν 78MaSi	ν 80Lee	ν 82NuGl	ν This work
0	—	65 647	—	65 693
1	—	66 520	—	66 691
2	—	67 376	—	67 434
3	—	68 227	68 217	68 194
4	69 085	69 085	69 080	68 997
5	69 950	69 940	69 935	69 864
6	70 750	70 746	70 746	70 701
7	71 515	71 515	71 515	71 516
8	72 280	72 290	72 296	72 276
9	73 035	73 041	73 046	73 047
10	73 770	73 746	73 741	73 738
11	74 570	74 555	74 566	74 600
12	75 370	75 432	75 438	75 501
13	—	76 214	76 220	76 294
14	—	—	76 938	76 966
RMSE	64	68	47

---

Fig. 7 compares a computed spectrum of the 140 nm band, calculated at 300 K with 0 < J ≤ 50, with the experiments. A partition function of 899.7 calculated using this rotational interval is in agreement within 99% with the values reported in the ExoMol database as taken from the Harris line list for HCN⁹¹ at this temperature, suggesting a near complete solution. The ExoMol partition function is based on rovibrational energies with J_max = 60. The magnitudes of the cross-sections are lower than those measured by 78MaSi and 80Lee suggesting that the 16ChRoAg transition dipoles are too small, while the cross-sections from 82NuGl are consistently higher than those in all experiments. For higher vibrational levels, we observe that the computed spectra show a higher level of rotational splitting than experiments. The vibrational progression of v₂ = 8 at 137.74 nm is spread along a 1 nm interval with a large splitting between the rotational components.

Fig. 7 Photoabsorption spectra of the 3 A′ ← 1 A′ transition. Our computed spectrum (0 < J < 25) is in black, experimental spectra are in blue for 78MaSi,⁸⁸ green for 80Lee,⁴⁵ and gold for 82NuGl.⁸⁹ The photodissociation region is found at wavelengths shorter than 138.84 nm.

The 2 A′ state corresponds to one component of the ¹Δ state at linear configurations, making transitions to this state forbidden at linear geometries and weak for bent structures. The intensities are of the order of 10⁻²⁰ cm⁻¹, being between 2 and 3 order of magnitudes lower than the other transitions. Fig. 8 shows our computed spectrum between 145 and 175 nm when using 1000 roots. A small vibrational progression is observed for wavelengths longer than 155 nm, followed by a bound-free region below this limit.

Fig. 8 Photoabsorption spectra of the 2A′ ← 1 A′ transition. The spectrum is computed at T = 0 K. The photodissociation regime is found for wavelengths below 155 nm.

Table 4 and Fig. 9 show our results for the 4 A′ ← 1 A′ and the 5 A′ ← 1 A′ bands compared to the experiments from 81NaKoOz and 82NuGl. Our peak assignments are made by performing calculations at T = 0 K. The experiments, however, report only the vibrational levels without information on the rotational levels. A small shift between the energy levels can be expected, as at room temperature the maximum of the cross-section will be found for higher rotational states. 81NaKoOz identifies 6 vibrational states within the two electronic state, and we found that our levels are consistently higher than theirs with an RMSE of 35.1 cm⁻¹. We get better agreement with the 17 levels measured by 82NuGl, with a RMSE of 14.0 cm⁻¹. A larger disagreement is found for the v₂ = 2 state, with a difference of 157 cm⁻¹. Excluding this state, the RMSE decreases to 10.0 cm⁻¹. However, we get poor agreement with the magnitudes of the cross-sections; our computed magnitudes are lower by a factor of 5–10 in the 100–130 nm region. It is likely that the computed TDMSs are too small.

Table 4 Vibrational band centers 4 A′ ← 1 A′ and 5 A′ ← 1 A′ bands: Comparison between theory (this work) and experiment, from 82NuGl⁸⁹ and 81NaKoOz⁹⁰

v 1 v 2 v 3	81NaKoOz	82NuGl	This work
4 A′ ← 1 A′
000	82 100	82 217	82 243
020	—	83 222	83 227
001	83 960	83 991	84 041
002	—	85 948	86 000
5 A′ ← 1 A′
000	88 980	89 063	89 013
010	—	89 718	89 666
020	90 170	90 400	90 243
001	90 930	90 992	91 025
011	—	91 659	91 706
100	92 070	92 039	92 113
021	—	92 396	92 508
110	—	92 653	92 682
002	—	92 963	92 936
031	—	93 101	93 090
120	—	93 396	93 384
130	—	94 171	94 159
200	—	95 057	95 057

---

Fig. 9 Bound–bound photoabsorption spectra of the 4 A′ ← 1 A′ (black) and 5 A′ ← 1 A′ (pink) transitions versus the experiments using a logarithmic y-scale. Experimental spectra are in green for 80Lee, orange for 81 NaKoOz and gold for 82NuGl. The computed spectra are consistently lower than that measured by the experiments.

3.2 Integrated intensities

The calculation of the integrated intensity (I) is a second method to assess the quality of our calculations. The method itself is the integration of the photodissociation cross-sections over a specific wavelength range (λ₁ and λ₂):

The MPI-Mainz UV/VIS Spectral Atlas (ref. 92)‡ contains 4 different measured electronic spectra of HCN. The intensities of the 3 A′ ← 1 A′ are integrated from the measurement from 78MaSi⁸⁸ and 80Lee.⁴⁵ Upper state excitations are measured from 80Lee, 81NaKoOz⁹⁰ and 82NuGl.⁸⁹ The Mainz database simply reports the overall spectrum over the experimental wavelength range, making it impossible to make a one to one comparison with our results in terms of individual transitions.

The integrated intensities for the 3 A′ ← 1 A′ bands from our model (1.57 × 10⁻¹⁴ cm molecule⁻¹) and the spectrum extracted from 16ChRoAg⁵⁹ (1.40 × 10⁻¹⁴ cm molecule⁻¹) are within 12%, showing overall agreement between the two computational models. For these calculations, only the J = 1 ← 0 transitions are considered. At room temperature (300 K), the calculated band intensity, with J varying from 0 to 25, is 3.01 × 10⁻¹⁴ cm molecule⁻¹. Experiments from 78MaSi, 80Lee, and 82NuGl show higher intensities than the computed ones, by a factor of 1.5 or 6, depending on the source. Values are reported in Table 5.

Table 5 Integrated intensities from experiments and calculations for HCN electronic transitions. Lee⁴⁵ and MacPherson & Simons⁸⁸

Source	T (K)	J	I (cm molecule^−1)
2 A′ ← 1 A′
Computed
This work	0	0,1	5.56 × 10^−17
This work	300	0,20	8.36 × 10^−18
16ChRoAg^59	—	0,1	4.17 × 10^−17
3 A′ ← 1 A′
Computed
This work	0	0,1	1.57 × 10^−14
This work	300	0,25	3.01 × 10^−14
16ChRoAg^59	—	0,1	1.40 × 10^−14
Experiments
80Lee^45	Room temperature	—	8.33 × 10^−14
82NuGl^89	Room temperature	—	1.86 × 10^−13
78MaSixx^88	Room temperature	—	4.41 × 10^−14
4 A′ ← 1 A′
Computed
This work	0	0,1	2.18 × 10^−14
This work	300	0,20	2.67 × 10^−14
5 A′ ← 1 A′
Computed
This work	0	0,1	6.07 × 10^−15
This work	300	0,20	8.71 × 10^−15
Upper states
Experiments
80Lee^45	Room temperature	—	2.50 × 10^−12
82NuGl^89	Room temperature	—	1.08 × 10^−12
81NaKoOz^90	Room temperature	—	8.20 × 10^−13

---

The integrated intensity of the 2 A′ ← 1 A′ band for our model agrees with the value from 16ChRoAg. The intensity decreases with the increase of temperature and rotational states, from 5.56 × 10⁻¹⁷ cm molecules at 0 K to 8.36 × 10⁻¹⁸ at 300 K.

There is no experimental measurement for individual higher states; rather, they are measured collectively. The intensities are of the order of magnitude of 10⁻¹² cm molecule⁻¹. These values are higher than our computed values by a factor between 30 and 90, depending on the experimental source. The cause of the underestimation of the cross-sections and intensities with respect to the experiment was ascribed by 16ChRoAg⁵⁹ to excited states possibly stealing intensity from dark states, not considered in their work or ours, through coupling to other states. It would seem more likely that these discrepancies arise from the computed PESs and TDMs: the electronic intensities are very sensitive to the position and shape of the upper state; coupled with a possible underestimation of the transition dipole moments.

The causes of the underestimation of the cross-sections and intensities with respect to the experiment was ascribed by 16ChRoAg for two possible reasons:⁵⁹ these excited states might steal intensity from dark states, not considered in this work, through coupling from other states. Other possible reasons of these discrepancies are in the computed PESs and TDMs: the electronic intensities are very sensitive to the position and shape of the upper state; coupled with a possible underestimation of the transition dipole moments.

3.3 Temperature dependence of the 3 A′ ← 1 A′ band

As in the case of diatomic molecules, the increase in temperature leads to a change in the shape of the cross-sections. In order to understand this effect, the cross-sections were computed at 300, 1000 and 2000 K and are plotted in Fig. 10 in linear and logarithmic forms. The calculated partition function of 50 681.8 at 2000 K is 99% of the value in the Harris line list for HCN⁹¹ from the ExoMol database. The increase of temperature corresponds to an increase of the background continuum, the well separated vibrational progression up to 2000 K. The maximum of the cross-section, corresponding to v₂ = 7 peak, is red shifted by 1.79 nm (890 cm⁻¹) from 300 K to 1000 K, and it remains constant up to 2000 K. The cross-section height decreases by a factor of 6 when the temperature reaches 2000 K from 300 K. These values are collected in Table 6. However, importantly for modelling applications, the threshold to photodissociation cross-section drops significantly with temperature, as observed previously in our study of diatomics,¹⁸ which can lead to very large increases in the effective rate of dissociation.

Fig. 10 Photoabsorption spectra of the 3 A′ ← 1 A′ transition at 300, 1000, and 2000 K. The main plot is in a linear scale, while the inset is the same plot but in the logarithmic scale.

Table 6 Temperature dependence of various spectroscopic properties (σ_max, λ_max) of the 3 A′ ← 1 A′ transition

T (K)	σ max (cm^2 molecule^−1)	λ max (nm)
300	1.67 × 10^−17	140.92
1000	8.22 × 10^−18	142.71
2000	3.04 × 10^−18	142.71

---

4 Conclusions

In this work, we explore the possibility of studying the direct photodissociation via dissociative potential surfaces and indirect photodissociation via quasi-bound excited states of triatomic molecules. The distinction between these two mechanisms is less clear cut in polyatomics than in diatomics. Our approach is based on solving the time-independent Schrödinger equation using variational nuclear motion calculations provided by adapting existing programs. We chose HCN as a test molecule, as it is a well-studied molecule experimentally and various computational studies. We performed the calculations using EVEREST, a program suite developed from DVR3D with some advantages when studying unbound states. This is because the latter code utilises a sinc-DVR basis, which provides better coverage of the dissociating coordinates compared to spherical oscillators, which tend to become increasingly sparse at large interatomic separations.

Our method is able to reproduce the experiments and the original work from 16ChRoAg.⁵⁹ The discrepancies between experiments and our calculations can be explained by our use of an adiabatic description of the PESs and issues with the calculated transition dipole moments. Our cross-sections and intensities, like as in the previous study by 16ChRoAg, are underestimated with respect to the experiments. In their paper, 16ChRoAg mention as possible reasons for these discrepancies the lack of couplings in the model and the difficulty of computing accurate transition dipole moments involving electronically excited states.

The temperature dependence of the cross-section has been tested on the 3 A′ ← 1 A′ transitions at three different temperatures. A shift in the cross-section peak and a decrease of the cross-section height plus a lower of the photodissociation threshold are observed with increasing temperatures, as found for diatomic molecules.^14,18

The next steps will be to produce a complete photodissociation model for HCN, performing electronic structure calculations for the states not included in this work and refining the available potential energy and transition dipole surfaces to recover the low-temperature experimental results. This model will then be used to provide cross-sections over an extended range of molecular temperatures. We plan to use our newly developed method to study other important molecules such as H₂O and H₂S. Non-adiabatic couplings between potential surfaces can provide an important mechanism for photodissociation. We note that EVEREST allows the implementation of different forms of surface couplings, including non-adiabatic effects. These features have been successfully tested on bound–bound electronic transitions of CaOH, which are characterised by both Renner–Teller and spin–orbit coupling effects, as demonstrated by Owens et al.⁷⁷ This is a supplementary reason for our choice of EVEREST. Calculations on the photodissociation of H₂S which include the explicit treatment of non-adiabatic couplings using EVEREST are currently in progress.

Data availability

The data for this article are contained in the article: numerical values for the figures are available from the authors. The documented codes, ExoSmooth and DVR3D, are available from https://github.com/exomol.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank professor Octavio Roncero for providing the potential energy surfaces and transition dipole moment surfaces. This work was funded by the ERC Advanced Investigator Project 883830 and the STFC Project ST/Y001508/1.

References

1. V. H. J. Clark and S. N. Yurchenko, Phys. Chem. Chem. Phys., 2021, 23, 11990–12004.
2. S.-M. Tsai, et al. , Nat. Astron., 2023, 617, 483–487.
3. O. Venot, M. Rocchetto, S. Carl, A. R. Hashim and L. Decin, Astrophys. J., 2016, 830, 77.
4. M. A. Badhan, E. T. Wolf, R. K. Kopparapu, G. Arney, E. M.-R. Kempton, D. Deming and S. D. Domagal-Goldman, Astrophys. J., 2019, 887, 34.
5. B. Fleury, M. S. Gudipati, B. L. Henderson and M. Swain, Astrophys. J., 2019, 871, 158.
6. N. K. Lewis, H. R. Wakeford, R. J. MacDonald, J. M. Goyal, D. K. Sing, J. Barstow, D. Powell, T. Kataria, I. Mishra, M. S. Marley, N. E. Batalha, J. I. Moses, P. Gao, T. J. Wilson, K. L. Chubb, T. Mikal-Evans, N. Nikolov, N. Pirzkal, J. J. Spake, K. B. Stevenson, J. Valenti and X. Zhang, Astrophys. J., 2020, 902, L19.
7. O. Venot, N. Fray, Y. Bénilan, M.-C. Gazeau, E. Hébrard, G. Larcher, M. Schwell, M. Dobrijevic and F. Selsis, Astron. Astrophys., 2013, 551, A131.
8. R. Valiev, A. Berezhnoy, I. Gritsenko, B. Merzlikin, V. Cherepanov, T. Kurten and C. Wöhler, Astron. Astrophys., 2020, 633, A39.
9. A. N. Heays, A. D. Bosman and E. F. van Dishoeck, Astron. Astrophys., 2017, 602, A105.
10. A. Noelle, A. C. Vandaele, J. Martin-Torres, C. Yuan, B. N. Rajasekhar, A. Fahr, G. K. Hartmann, D. Lary, Y.-P. Lee, P. Limao-Vieira, R. Locht, K. McNeill, J. J. Orlando, F. Salama and R. P. Wayne, J. Quant. Spectrosc. Radiat. Transfer., 2020, 253, 107056.
11. S. Y. Grebenshchikov, J. CO2 Util., 2016, 15, 32–40.
12. J. Tennyson, S. N. Yurchenko, A. F. Al-Refaie, E. J. Barton, K. L. Chubb, P. A. Coles, S. Diamantopoulou, M. N. Gorman, C. Hill, A. Z. Lam, L. Lodi, L. K. McKemmish, Y. Na, A. Owens, O. L. Polyansky, T. Rivlin, C. Sousa-Silva, D. S. Underwood, A. Yachmenev and E. Zak, J. Mol. Spectrosc., 2016, 327, 73–94.
13. J. Tennyson, S. N. Yurchenko, A. F. Al-Refaie, V. H. J. Clark, K. L. Chubb, E. K. Conway, A. Dewan, M. N. Gorman, C. Hill, A. E. Lynas-Gray, T. Mellor, L. K. McKemmish, A. Owens, O. L. Polyansky, M. Semenov, W. Somogyi, G. Tinetti, A. Upadhyay, I. Waldmann, Y. Wang, S. Wright and O. P. Yurchenko, J. Quant. Spectrosc. Radiat. Transfer, 2020, 255, 107228.
14. M. Pezzella, S. N. Yurchenko and J. Tennyson, Phys. Chem. Chem. Phys., 2021, 23, 16390–16400.
15. S. N. Yurchenko, L. Lodi, J. Tennyson and A. V. Stolyarov, Comput. Phys. Commun., 2016, 202, 262–275.
16. I. I. Mizus, L. Lodi, J. Tennyson, N. F. Zobov and O. L. Polyansky, J. Mol. Spectrosc., 2022, 368, 111621.
17. S. N. Yurchenko, A. F. Al-Refaie and J. Tennyson, Astron. Astrophys., 2018, 614, A131.
18. M. Pezzella, S. N. Yurchenko and J. Tennyson, Mon. Not. R. Astron. Soc., 2022, 514, 4413–4425.
19. Z. Qin, T. Bai and L. Liu, Mon. Not. R. Astron. Soc., 2022, 516, 550–561.
20. L. E. Snyder and D. Buhl, Astrophys. J., 1971, 163, L47.
21. T. Hirota, S. Yamamoto, H. Mikami and M. Ohishi, Astrophys. J., 1998, 503, 717–728.
22. J. T. van Loon, J. R. Marshall, M. Cohen, M. Matsuura, P. R. Wood, I. Yamamura and A. A. Zijlstra, Astron. Astrophys., 2006, 447, 971–989.
23. J. S. Carr and J. R. Najita, Astrophys. J., 2011, 733, 102.
24. S. Bruderer, D. Harsono and E. F. van Dishoeck, Astron. Astrophys., 2015, 575, A94.
25. T. Hirota, S. Yamamoto, K. Kawaguchi, A. Sakamoto and N. Ukita, Astrophys. J., 1999, 520, 895–900.
26. M. Lippi, G. L. Villanueva, M. A. DiSanti, H. Boehnhardt, M. J. Mumma, B. P. Bonev and D. Prialnik, Astron. Astrophys., 2013, 551, A51.
27. E. S. Wirström, M. S. Lerner, P. Källström, A. Levinsson, A. Olivefors and E. Tegehall, Astron. Astrophys., 2016, 588, A72.
28. E. Lellouch, P. N. Romani and J. Rosenqvist, Icarus, 1994, 108, 112–136.
29. E. Lellouch, B. Bezard, D. F. Strobel, G. L. Bjoraker, F. M. Flasar and P. N. Romani, Icarus, 2006, 184, 478–497.
30. E. Lellouch, M. Gurwell, B. Butler, T. Fouchet, P. Lavvas, D. F. Strobel, B. Sicardy, A. Moullet, R. Moreno, D. Bockelee-Morvan, N. Biver, L. Young, D. Lis, J. Stansberry, A. Stern, H. Weaver, E. Young, X. Zhu and J. Boissier, Icarus, 2017, 286, 289–307.
31. A. Tsiaras, M. Rocchetto, I. P. Waldmann, G. Tinetti, R. Varley, G. Morello, E. J. Barton, S. N. Yurchenko and J. Tennyson, Astrophys. J., 2016, 820, 99.
32. E. M. Molter, C. A. Nixon, M. A. Cordiner, J. Serigano, P. G. J. Irwin, N. A. Teanby, S. B. Charnley and J. E. Lindberg, Astrophys. J., 2016, 152, 42.
33. L. Rezac, A. A. Kutepov, A. Faure, P. Hartogh and A. G. Feofilov, Astron. Astrophys., 2013, 555, A122.
34. D. Dubois, N. Carrasco, M. Petrucciani, L. Vettier, S. Tigrine and P. Pernot, Icarus, 2019, 317, 182–196.
35. S. E. van Terwisga, E. F. van Dishoeck, P. Cazzoletti, S. Facchini, L. Trapman, J. P. Williams, C. F. Manara, A. Miotello, N. van der Marel, M. Ansdell, M. R. Hogerheijde, M. Tazzari and L. Testi, Astron. Astrophys., 2019, 623, A150.
36. F. Long, A. D. Bosman, P. Cazzoletti, E. F. van Dishoeck, K. I. Oberg, S. Facchini, M. Tazzari, V. V. Guzman and L. Testi, Astron. Astrophys., 2021, 647, A118.
37. M. Jin, J.-E. Lee and K.-T. Kim, Astrophys. J., Suppl. Ser., 2015, 219, 2.
38. J. Bublitz, J. H. Kastner, M. Santander-García, V. Bujarrabal, J. Alcolea and R. Montez, Astron. Astrophys., 2019, 625, A101.
39. J. Bublitz, J. H. Kastner, P. Hily-Blant, T. Forveille, M. Santander-García, J. Alcolea and V. Bujarrabal, Astron. Astrophys., 2022, 659, A197.
40. D. D. Davis and H. Okabe, J. Chem. Phys., 1968, 49, 5526–5531.
41. A. Mele and H. Okabe, J. Chem. Phys., 1968, 51, 4798–4808.
42. A. Meenakshi and K. K. Innes, J. Chem. Phys., 1986, 84, 6550–6554.
43. J. Guo, R. Eng, T. Carrington and S. V. Filseth, J. Chem. Phys., 2000, 112, 8904–8909.
44. S. Tatematsu and K. Kuchitsu, Bull. Chem. Soc. Jpn., 1977, 50, 2896–2898.
45. L. C. Lee, J. Chem. Phys., 1980, 72, 6414–6421.
46. C. R. Bucher and K. K. Lehmann, Chem. Phys. Lett., 1998, 294, 173–180.
47. U. G. Jørgensen, J. Almlöf, B. Gustafsson, M. Larsson and P. Siegbahn, J. Chem. Phys., 1985, 83, 3034–3042.
48. T. van Mourik, G. J. Harris, O. L. Polyansky, J. Tennyson, A. G. Császár and P. J. Knowles, J. Chem. Phys., 2001, 115, 3706–3718.
49. J. M. Bowman, S. Irle, K. Morokuma and A. Wodtke, J. Chem. Phys., 2001, 114, 7923–7934.
50. A. J. C. Varandas and S. P. J. Rodrigues, J. Phys. Chem. A, 2006, 110, 485–493.
51. R. Dawes, A. F. Wagner and D. L. Thompson, J. Phys. Chem. A, 2009, 113, 4709–4721.
52. V. Y. Makhnev, A. A. Kyuberis, O. L. Polyansky, I. I. Mizus, J. Tennyson and N. F. Zobov, J. Mol. Spectrosc., 2018, 353, 40–53.
53. M. Perić, H. Dohmann, S. Peyerimhoff and R. J. Buenker, ZPhys-e.D: At., Mol. Clusters, 1987, 5, 65–75.
54. D. Xu, D. Xie and H. Guo, Chem. Phys. Lett., 2001, 345, 517–524.
55. D. Xu, D. Xie and H. Guo, J. Chem. Phys., 2002, 116, 10626–10635.
56. D. Xu, H. Guo and D. Xie, J. Theor. Comput. Chem., 2003, 02, 639–648.
57. D. Xu, D. Xie and H. Guo, J. Phys. Chem. A, 2002, 106, 10174–10183.
58. M. K. Nayak, R. K. Chaudhuri and S. N. L. G. Krishnamachari, J. Chem. Phys., 2005, 122, 184323.
59. A. Chenel, O. Roncero, A. Aguado, M. Agúndez and J. Cernicharo, J. Chem. Phys., 2016, 144, 144306.
60. A. Aguado, O. Roncero, A. Zanchet, M. Agúndez and J. Cernicharo, Astrophys. J., 2017, 838, 33.
61. M. S. Priyadarshini, S. M. Jo, S. Venturi, D. W. Schwenke, R. L. Jaffe and M. Panesi, J. Phys. Chem. A, 2022, 126, 8249–8265.
62. H.-D. Meyer, U. Manthe and L. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
63. M. H. Beck, A. Jäckle, G. A. Worth and H. D. Meyer, Phys. Rep., 2000, 324, 1–105.
64. G. A. Worth, H.-D. Meyer, H. Köppel, L. S. Cederbaum and I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569–606.
65. U. Manthe, H. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213.
66. U. Manthe, H. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 9062–9071.
67. U. Manthe and A. D. Hammerich, Chem. Phys. Lett., 1993, 211, 7–14.
68. A. D. Hammerich, U. Manthe, R. Kosloff, H. Meyer and L. S. Cederbaum, J. Chem. Phys., 1994, 101, 5623–5646.
69. S. Han, C. Xie, X. Hu, D. R. Yarkony, H. Guo and D. Xie, J. Phys. Chem. Lett., 2023, 14, 10517–10530.
70. R. Schinke, Photodissociation Dynamics, Cambridge University Press, 1993.
71. G. B. Mitev, S. N. Yurchenko and J. Tennyson, J. Chem. Phys., 2024, 160, 144110.
72. J. Tennyson, M. Pezzella, J. Zhang and S. N. Yurchenko, RAS Tech. Instr., 2023, 2, 231–237.
73. S. N. Yurchenko, W. Szajna, R. Hakalla, M. Semenov, A. Sokolov, J. Tennyson, Y. Pavlenko and M. R. Schmidt, Mon. Not. R. Astron. Soc., 2024, 527, 9736–9756.
74. J. Tennyson, M. A. Kostin, P. Barletta, G. J. Harris, O. L. Polyansky, J. Ramanlal and N. F. Zobov, Comput. Phys. Commun., 2004, 163, 85–116.
75. A. O. Mitrushchenkov, J. Chem. Phys., 2012, 136, 024108.
76. E. J. Zak and J. Tennyson, J. Chem. Phys., 2017, 147, 094305.
77. A. Owens, A. Mitrushchenkov, S. N. Yurchenko and J. Tennyson, Mon. Not. R. Astron. Soc., 2022, 516, 3995–4002.
78. D. T. Colbert and W. H. Miller, J. Chem. Phys., 1992, 96, 1982–1991.
79. G. C. Groenenboom and D. T. Colbert, J. Chem. Phys., 1993, 99, 9681–9696.
80. J. Tennyson and B. T. Sutcliffe, J. Chem. Phys., 1982, 77, 4061–4072.
81. J. Tennyson and B. T. Sutcliffe, J. Mol. Spectrosc., 1983, 101, 71–82.
82. G. B. Mitev, R. P. Brady, O. Smola, J. Tennyson and S. N. Yurchenko, Phys. Chem. Chem. Phys., 2024 In preparation.
83. J. Tennyson and B. T. Sutcliffe, Mol. Phys., 1986, 58, 1067–1085.
84. J. Tennyson, C. Hill and S. N. Yurchenko, 6^th international conference on atomic and molecular data and their applications ICAMDATA-2012, 2013, pp. 186–195.
85. A. U. Hazi and H. S. Taylor, Phys. Rev. A: At., Mol., Opt. Phys., 1970, 1, 1109–1120.
86. Z. Bacic and J. Simons, J. Phys. Chem., 1982, 86, 1192–1200.
87. V. A. Mandelshtam, T. R. Ravuri and H. S. Taylor, Phys. Rev. Lett., 1993, 70, 1932–1935.
88. M. T. Macpherson and J. P. Simons, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 1965–1977.
89. J. A. Nuth and S. Glicker, J. Quant. Spectrosc. Radiat. Transfer, 1982, 28, 223–231.
90. T. Nagata, T. Kondow, Y. Ozaki and K. Kuchitsu, Chem. Phys., 1981, 57, 45–53.
91. R. J. Barber, J. K. Strange, C. Hill, O. L. Polyansky, G. C. Mellau, S. N. Yurchenko and J. Tennyson, Mon. Not. R. Astron. Soc., 2014, 437, 1828–1835.
92. H. Keller-Rudek, G. K. Moortgat, R. Sander and R. Sörensen, Earth Syst. Sci. Data, 2013, 5, 365–373.

---

Footnotes

† Present address: Dipartimento di Fisica e Geologia, Universitá di Perugia, 06123 Perugia, Italy.

‡ https://www.uv-vis-spectral-atlas-mainz.org/uvvis

---