A method for calculating temperature-dependent photodissociation cross sections and rates

The destruction of molecules by photodissociation plays a major role in many radiation-rich environments, including the evolution of the atmospheres of exoplanets, which often exist close to UV-rich stars. Most current photodissociation calculations and databases assume T = 0 K, which is inadequate for hot exoplanets and stars. A method is developed for computing photodissociation spectra of diatomic molecules as a function of temperature exploiting bound state variational nuclear motion program Duo and post-processing program ExoCross. Discrete transition intensities are spread out to represent a continuous photodissociation spectrum either by Gaussian smoothing or by averaging calculations over a range of different grid sizes. Our approach is tested on four different chemical species (HCl, HF, NaCl and BeH), showing its ability to reproduce photodissociation cross sections and rates computed with other approaches and experiment. The temperature dependence of photodissociation cross sections and rates is studies showing strong temperature variation of the photodissociation cross sections.


Introduction
Photochemistry substantially impacts the atmospheric composition of planets and exoplanets with consequences for the chemical compositions, radiative transfer, thermal structure, and dynamics of the atmospheres. This is particularly true for the many exoplanets that have been discovered orbiting near their host stars as these planets exists in UV-rich environments. [1][2][3][4] Of course, it is exactly the planets which experience high UV fluxes where the molecules are also hot and hence vibrationally and rotationally excited. Modelling and understanding the atmospheres of such planets therefore requires temperature-dependent photodissociation cross sections. Measurements of photodissociation cross sections of molecules at higher temperatures have been performed 1 but these studies struggle to reach the temperatures needed for the top of hot atmospheres (T 4 1000 K).
Current state-of-the-art for calculation of photodissociation cross sections for astronomical studies often use simplified (harmonic) ground state wavefunctions. 5 This model is appropriate for the cold molecules such as those found in the interstellar medium (ISM) but inadequate for hot environments such as the atmospheres of exoplanets. Here we present a novel methodology aimed at resolving this problem. Similarly, there are standard databases of photodissociation cross sections studies of the ISM 6,7 but these only contain data for molecules at interstellar temperatures, often assumed to be 0 K.
There is a long history of theoretical treatments of photodissociation, 8 but it is only recently that cross section calculations have begun to seriously consider the effects of temperature 9 and even then the effects of rotational excitation appears to have been largely ignored.
The ExoMol project was designed to produce comprehensive line lists of bound-bound transitions for molecules in hot atmospheres. 10 The ExoMol database provides such line lists for a large range of molecules deemed to be important in exoplanets and elsewhere. 11 As part of the ExoMol project a series of nuclear motion codes have been developed or enhanced to give results which are both comprehensive and accurate. 12 Here we concentrate on one of the programs, Duo. 13 Duo solves the bound-bound diatomic nuclear motion problem by explicit solution of the nuclear motion Schrödinger equation and allows for treatment of spin-orbit and other coupling effects 14 as well avoided and allowed curve crossings. 15 The treatment we propose here is based on extending Duo to treat the bound-free problem of photodissociation. This will allow temperature effects to be fully captured; our procedure can also provide the requisite data for modelling the effect of photodissociation in environments where non-local thermodynamic equilibrium (non-LTE) is important. Such data are not obtainable with current experimental procedures. We note that photodissociation itself is likely to prove to be a major driver of non-LTE regimes.
Duo has already been adapted for the study of continuum (free-free) states within an R-matrix formalism with a particular focus on the treatment of ultracold collisions. 16 Here we propose a rather more radical approach where the continuum is only modelled using a (finite) inner region and photodissociation cross sections or rates are extracted from the results. This approach avoids the need for computationally expensive treatment of the long-range wavefunctions making studies over many rotational states easy and fast. The approach allows the full photoabsorption problem, i.e., bound-bound and bound-free transitions, to be treated within a single formalism and on an equal footing.

Theoretical background
Our treatment of the diatomic photoabsorption problem, which includes both bound-bound transitions and bound-free (photodissociation) transitions is based on the diatomic code Duo. 13 Duo uses a variational procedure to find solutions to the multistate rovibronic nuclear motion problem allowing for treatment of spin-orbit and other couplings. 14 Duo has been extensively used to provide accurate line lists for challenging bound-bound problems. [17][18][19][20][21][22] The present implementation uses Duo to solve the nuclear motion Schrödinger equation and then post-processing program ExoCross 23 to produce photodissociation cross sections and rates, and, optionally, to distinguish between bound-bound and bound-free transitions. Duo initially uses a grid basis set to solve the rotationless one-dimensional Schrödinger equation separately for each electronic state G: producing a set of vibrational wavefunctions f G n . Here V G (r) is the corresponding potential energy curve (PEC), m is the reduced mass of the molecule, r the interatomic distance, and n the vibrational quantum number. In the present implementation solutions can be obtained using a grid based on a sinc DVR (discrete variable representation) 24 or on Lobatto shape functions 16,25 or a fivepoint finite differences to derive the kinetic energy operator.
The transition dipole moment curves, T G 0 ;G 00 n 0 ;n 00 ðrÞ between two vibronic states f G 0 n 0 and f G 00 n 00 are expressed as: with m G 0 ;G 00 ðrÞ as the electronic dipole moment vector and are used for evaluating the intensity of the transitions.
The vibrational wavefunction are then combined with suitable angular functions in a Hund's case a representation, to provide a basis set for each total angular momentum quantum number, J. These basis sets are used to solve the full Hamiltonian, that couples rovibronic states c J i belonging to different electronic states and different values of the angular momenta J. Additional terms can be added to take into account non-adiabatic couplings between curves, 15 and allow for spinorbit and other similar couplings. Only transitions f ' i that obey the electric dipole moment selection rules parity changes and are allowed.
The intensity I f'i of a given bound-bound transition also depends on the rotational number of the initial ( J i ) and final ( J f ) states and on the temperature (T): where A fi is the Einstein-A coefficient (s À1 ) computed using the Duo rovibronic wavefunctions |c J i i,ñ fi is the transition wave- is the second radiations constant (cm K À1 ); h is the Planck constant, c the speed of light, k B the Boltzmann constant; g tot i is the total nuclear statistical weight factor where g ns i is nuclear spin statistical weight; Q(T) is the partition function defined as a sum over bound states The intensity is the integral of the cross section s fi over an absorption line: By introducing a line profile fñ fi (ñ), s fi (ñ) can be defined as s fi (ñ) = I f'i fñ fi (ñ), where fñ fi (ñ) is an integrable function which is normalized to unity. The bound-free photodissociation process is characterized by the excitation from a bound electronic state, usually the electronic ground state, to an unbound rovibronic level of an excited state. The radial wavefunctions of these dissociative states are described by a sinusoidal wavefunction at the asymptotic limit. Here we adapt eqn (4) to cover both bound-bound and bound-free processes. We note that the self-absorption term in eqn (4), given by Àe Àc2ñfi/T , is probably not needed for bound-free transitions 26 but in practice will be negligible for the short wavelength processes considered here. A future refinement will be to remove this term for bound-free processes.

Continuum cross section calculations
Duo is designed to provide discrete solutions for bound electronic systems of diatomics. Here we present a robust approach to use Duo for computing temperature-dependent photo-dissociation spectra of diatomics representing bound-free transitions. To this end, a coupled set of Schrödinger equations for the system containing bound and unbound PECs is solved on the basis of bound vibrational functions f G from eqn (1). The discrete eigenvalues Ẽ i and eigenfunctions c J i are then used to generate line intensities via eqn (4). This gives a photodissociation spectrum which is represented by clusters of discrete, high intensity lines, separated by regions where there is no intensity, see Fig. 1. To recover the continuum nature of the spectrum we apply a smoothing function to the cross sections, computed using the ExoCross program. 23 We have tested two smoothing functions included into the SciPy 27 package. The first method consists in interpolating the spectrum with knots equally distributed along the wavelengths. The second method consist of applying a normalized Gaussian smoothing function to each grid point. The Gaussian line profile is given by: 28 whereñ fi is the line centre and a D is the Gaussian half width at half maximum (HWHM). The HWHM of the Gaussian line profile is used for regulating the cross section height. Optimal values for the HWHMs depend strongly on the molecule under analysis, with values between 2000 cm À1 and 3000 cm À1 for hydrogen halides (see Section 3.1), and 75 cm À1 in case of high temperature cross sections for BeH + (Section 3.3). The appropriate HWHM depends on several factors including the size and density of the grid used, which determines the number of discrete transitions given by Duo, and the temperature as the large number of lines at elevated temperatures leads to dense spectra making further smoothing beyond the initial 10 cm À1 broadening unnecessary. Both methods are designed to conserve the integrated intensity and the photodissociation rates with respect to the initial transition intensities.
Some excited potentials support a few bound vibrational states. These states manifest themselves in a photoabsorption spectrum as a series of bound-bound transitions at longer wavelengths than the transitions to those states which are responsible for photodissociation. It is therefore important to be able distinguishing bound-free transitions from boundbound ones. Here we adopt an approach that has some similarities to the stabilization method of Taylor and co-workers. [29][30][31] In our case, the Duo calculations are repeated using different grid sizes which are varied by a few tenths of an Å. Each calculation results in temperature-depended cross sections, which are then averaged to produce our final photodissociation cross sections. The individual cross sections are obtained using a Gaussian line profile of HWHM = 10 cm À1 . These repeated calculations smear out the bound-free transitions, possibly leading to results that are more easily smoothed to give continuum cross sections, but leave bound-bound transitions in the same place allowing them to be readily identified. Photodisocciation cross sections generated using the stabilization method are illustrated in Fig. 1.
If one only wants photodissociation cross sections or rates, bound-bound contributions need to be identified and discounted. Identification of bound-bound transitions is facilitated by using the stabilization method as they always occur with the same transition frequency when the box size is varied. The photodissociation spectrum can be recovered by calculating the overall photoabsorption spectrum and then subtracting the bound-bound transitions contribution. An alternative approach consists of summing photodissociation cross sections evaluated for each single state in turn excluding the bound contributions. Results of these two methods are compared in Section 3.3.

Photodissociation rates
For many purposes photodissociation rates are used instead of cross sections; for example, rates are used for modelling the abundance and the evolution of species in space. 32,33 The rates provide a useful quantity to test the validity of our approach for the molecules we study. The photodissociation rate k of a molecule dissociated by a field with a flux F(ñ) between the wavelengthsñ 1 andñ 2 is expressed as: There are several standard fluxes used to produce appropriate rates. Here we concentrate on the flux appropriate for the interstellar medium (ISM) since this is widely used by the databases to which we want to compare. Future work will consider a variety of stellar fluxes. Photodissociation rates are calculated using the interstellar radiation field (ISRF). [34][35][36] The ISRF has been fitted to an analytical expression for wavelengths between 91.2 nm and 200 nm, and was expressed as: where l is the wavelength in nm and it was later extended to 2000 nm using the expression: 3 Results Our approach is tested for three different system types. The first consists of the A 1 P ' X 1 S + photodissociation from the Fig. 1 Cross sections generated using a Gaussian with HWHM of 10 cm À1 for HCl and HF from a single run in Duo (black) and with the stabilization method (red). Our calculations are performed at 100 K.
vibrational ground state for two hydrogen halides: HCl and HF.
The second system is NaCl and its photodissociation as a function of temperature. With the third system, BeH + , we compare the results from our method with recently published calculations. All calculations are carried using a 8 core Intel local machine with 32 GB of RAM using available potential energy and transition dipole curves.

HCl and HF
HCl and HF are characterized by the presence of a repulsive A 1 P excited state. The A 1 P' X 1 S + electronic transition leads to the immediate photodissociation into two neutral fragments H( 2 S) + X( 2 P), where X is F or Cl. Photodissociation arising from these electronic transitions was chosen as an initial test of our methodology. For HCl we used the X 1 S + and the A 1 P potentials taken from Alexander et al. 37 and the A 1 P ' X 1 S + transition dipole moment from Givertz and Balint-Kurti. 38 For HF the X 1 S + , A 1 P potentials and the A 1 P ' X 1 S + transition dipole moments of HF are taken from Brown and Balint-Kurti. 39 For both molecules, vibrational wavefunctions were built for J = 0, between 0.5 Å and 3.0 Å. Transitions from the vibrational ground state of the X 1 S + state to the A 1 P state were considered for a temperature of T = 100 K, which required J up to 16 for HCl and up to 11 for HF. Calculations were also performed with the time-dependent Schrödinger code PHOTO, developed by Balint-Kurti et al., 40 which only considers states with J = 0 and is hence useful for a low temperature comparison. Fig. 1 presents photodissociation cross sections of HCl and HF computed using the stabilization approach for T = 100 K. The black lines in Fig. 1 show the intensity of the discretised transitions to the continuum obtained running a single calculation with Duo and ExoCross on a grid of 2001 points ranging from 0.5 Å to 3.00 Å. These spectra of HCl and HF consist in clusters of discrete lines, characterized by high values of s(ñ), separated by regions where s(ñ) = 0. The red curves in Fig. 1 show the spectra calculated with the stabilization method extending the original grid from 2.50 Å to 2.60 Å with steps of 0.001 Å, obtained by averaging 100 individual cross sections, each computed using the Gaussian line profile with HWHM = 10 cm À1 . The final spectrum consists of a discrete spectrum overlapped to a continuum background. The magnitude of the cross section for the discrete spectrum obtained with this procedure is 20 times smaller the peaks in the black curves of Fig. 1. Now we apply the Gaussian smoothing method to produce the T = 100 K photodissociation cross sections of HCl as described above. These are compared to experimental results [41][42][43] and to the results from PHOTO 40 with the numerical results are reported in Table 2 and plotted in Fig. 2. The Gaussian smoothing model with HWHM of 1800 cm À1 produces cross sections with s max = 3.55 Â 10 À18 cm 2 molecule À1 , while the interpolation scheme (with knots at every 5 nm) leads to a higher value of s max = 4.13 Â 10 À18 cm 2 molecule À1 , but within the upper limits reported by Inn. 41 The peak position (l max ) is overestimated by 3 Å with respect to PHOTO and experiments. Our photodissociation rates calculated for ISRF before and after smoothing give the same value of k = 2.29 Â 10 À10 s À1 at T = 100 K. The photodissociation rates differ between our model and PHOTO only by the 4%.
The experimental HCl photodissociation cross section shows an asymmetry at short wavelengths, see Fig. 2, due to a non-adiabatic coupling between the A 1 P and the C 1 P states. 43 The maximum of the cross section is found at l max = 153.0 nm, with s max = 3.53 Â 10 À18 cm 2 molecule À1 . The computer model of van Dishoeck et al. 44 overestimates l max by 2.3 nm with respect to the experiments and estimates a photodissociation rate of k = 2.1 Â 10 À10 s À1 . The time dependent Schrödinger code PHOTO estimates the photodissociation rate of k = 2.188 Â 10 À10 s À1 .
Three different types of basis functions were tested in Duo, corresponding to Sinc DVR, Lobatto and the 5 point finite differences; Table 1 shows the corresponding values of l max , s max as well as of the integrated intensity I. The calculations were performed on a grid of 601 points. The Lobatto wavefunctions lead to the largest difference from the other methods, especially in case of HF. For this molecule we observe a shift of the maximum cross section peak l max of 2.08 nm and an overestimation of the cross section maximum of the order of 18%. All further calculations in this work are performed using the sync DVR basis function.
There is a weak dependence of the cross section on the number of grid points. Increasing the size from 251 to 4001 points, for both HCl and HF, there is an increase of 0.6% in s max and of 0.3% in l max .
The theoretical and experimental photodissociation cross sections of HF are compared in Fig. 3. The shape of the HF photodissociation cross section is symmetric, with measured values of l max varying between 119.8 and 121.7 nm. 45,46 s max also varies greatly between the two experiments, from 6 Â 10 À18 cm 2 molecule À1 given by Hitchcock et al. 45 44 and the experimental data in cyan. 43 Our calculations are performed at 100 K.
Our l max for HF at T = 100 K agrees with Hitchcock et al. 45 The spectrum from Duo and ExoCross show a high energy tail at 90 nm (see Fig. 3), that is reproduced in the Gaussian smoothing model with a HWHM of 2700 cm À1 . Our photodissociation rates (1.173 Â 10 À10 s À1 ) agree within an uncertainty of 4% with the rates calculated by Brown and Balint-Kurti 39 and our calculations using PHOTO. Table 3 compares calculated and experimental values of l max , s max and k.
The temperature dependence of s max and k is reported in Fig. 4. For temperatures below 1000 K, s max is constant and it starts decreasing at higher temperatures. In the case of HCl, at 0.01 K s max = 3.53 Â 10 À18 cm 2 molecule À1 and decreases to s max = 4.77 Â 10 À19 cm 2 molecule À1 at 20 000 K for the Gaussian smoothing model. The same trend is observed for HF, with s max = 3.13 Â 10 À18 at 0.01 K and s max = 5.50 Â 10 À19 at 20 000 K. The bottom part of Fig. 4 shows that photodissociation rates for HCl and HF with the expected decrease when temperature increases. The figure shows near perfect agreement, within the 0.2%, between the rates calculated using the raw, unsmoothed data from ExoCross and the two smoothing models.
The temperature increase leads to an increase of population of the excited (bound) rovibrational states of the lower, ground electronic state. As a consequence, the photodissociation spectrum has a more dense character forming featureless spectral background even without extra smoothing applied. Fig. 5 shows the unsmoothed and Gaussian smoothed cross section for HCl at different temperatures: the two of them present the same shape for temperatures equal to or higher than 2000 K.

NaCl
Experimental photodissociation cross sections for NaCl have been reported for T = 300 K Na + + Cl À dissociation products, 47 Table 1 Comparison of results obtained using three different basis functions, corresponding to the sinc DVR, Lobatto wavefunctions, and five-point finite differences methods, for HCl and HF. Calculations are performed on 601 point grid at 100 K for an interatomic distance between 0.5 and 3 Å. The properties considered are the wavelength of the cross section peak (l max , nm), the cross section maximum (s max , cm 2 molecule À1 ), the integrated intensity (I, cm molecule À1 ), and the partition function (Q) at T= 100 K  Fig. 3 Cross sections for HF for computational models presented in Table 3. Gaussian smoothing results are in red, results of interpolation in green, and calculation performed with PHOTO in blue; the calculations of Brown and Balint-Kurti 39 are in violet. The Gaussian smoothing model is sensitive to the high energy tail shown in Fig. 1.   and at T = 1123 K for the neutral dissociation products. 48 The photodissociation spectrum above 200 nm comprises two contributions from the A 1 P ' X 1 S + and the B 1 S + 'X 1 S + transitions. At low temperatures, both bands form distinct, observable structures, as shown in Silver et al., 47 while they merge at higher temperatures: the cross section at 1123 K shows a maximum at l = 236 nm with s = 3.5 AE 0.3 Â 10 À17 cm 2 . 48 Our calculations are performed using the potential energy curves, dipole moments and couplings from the ExoMol study of Barton et al.,49 where they were used to calculate the X 1 S + state rovibrational spectrum. As part of the current work, we have produced components required for modelling the electronic A-X and B-X spectra using the MRCI/aug-cc-pVQZ ab initio level of theory as implemented in MOLPRO. 50 This includes the PECs X 1 S + , A 1 P and B 1 S + , the transition dipole moment curves A-X and B-X and an electronic angular momentum coupling curve between the A and B state. The active space was selected to be (10, 4, 4, 0) with (4, 2, 2, 0) closed orbitals.
For temperatures below 1000 K, transitions up to J = 100 are used; this is extended to J = 291 for higher temperatures. A 1001 point grid with points between 1.5 Å and 5 Å is used, selecting 60 vibrations states for the X 1 S + states, 150 for the A 1 P and 120 for the B 1 S + . The electronic angular momentum coupling curve L x is considered between the X 1 S + and the A 1 P states. The calculated photodissociation cross sections of NaCl for different temperatures are plotted in Fig. 6 with the corresponding l max and s max values tabulated in Table 4. A direct comparison between our calculations and experimental cross sections of NaCl at 1123 K are shown in Fig. 7. Our calculations reproduce the two peak structure in the NaCl cross sections and the disappearance of the lower energy peak at high temperatures. The photodissociation band A 1 P ' X 1 S + is very distinct at T = 100 K, T = 300 K, T = 500 K, turning a into shoulder at T = 750 K. It is submerged by the B 1 S + ' X 1 S + band at higher temperatures. For the analysis presented in Table 4, their contributions for T Z 1200 K were separated as indicated by an asterisk. The B 1 S + ' X 1 S + transition shows a blue shift of 1.71 nm and a decrease of the cross section from 19.1 Â 10 À16 to 2.67 Â 10 À17 over the 100-1500 K interval. A direct comparison between the cross sections at 300 K of Silver et al. 47 is not possible, due to the different final states. The photodissociation rate for the ISRF field shows almost no temperature dependence, passing from 9.24 Â 10 À10 s À1 at T = 100 K to 9.26 Â 10 À10 s À1 at T = 1500 K.
The experimental T = 1123 K cross section by Davidovits and Brodhead 48 (see Fig. 7) shows a feature from the C 1 P ' X 1 S + band at shorter wavelengths, which is not present in our model. The smoothed curve is characterized by l max = 233.82 nm versus the experimental value of l max = 236 nm with a cross section of 3.17 Â 10 À17 cm 2 molecule À1 , within the uncertainty range of the experimental results.

BeH +
Beryllium is the lightest stable nuclide not synthesized in the Big Bang, and it is a probe to study the early Universe structure and evolution. 51,52 The hydrogenated species BeH + finds application in plasma and nuclear physics. 53,54 The potential energy curves, transition dipole moments, spectroscopic data and photodissociation cross sections of BeH + covering the X 1 S + , A 1 S + , B 1 P, C 1 S + states at 1800 K, 4500 K, 10 000 K and 20 000 K have been reported by Xu et al. 55 and Yang et al. 56 The photodissociation spectrum of BeH + comprises two peaks, the first from the A 1 S + , B 1 P ' X 1 S + bands that correlates with the Be + (2p) + H(1s), and the second one from the C 1 S + ' X 1 S + band, which correlates with the Be(2s 2 ) + H + asymptote. All the three excited electronic states support bound vibrational states: 23 for A 1 S + , 13 for B 1 P, and 7 for C 1 S + . Fig. 6 NaCl cross sections calculated at different temperatures, T = 100 K black, T = 300 K red, T = 500 K green, T = 750 K blue, T = 1200 K orange and T = 1500 K violet. The A 1 P ' X 1 S + contribution and the vibrational structures disappear with the increase of temperature. Table 4 NaCl cross section and rates at different temperatures. Contributions deriving from the B1S + ' X 1 S + and A 1 P ' X 1 S + transitions, are separated. For T 4 1200 K, the low energy cross section is submerged by the high energy one; these cases are evidenced by an asterisk. The temperature is in Kelvin, l B'X max and l A'X max are in nm, s B'X max and s A'X max are cm 2 molecule À1 , and k in s À1 The coupled system of Schrödinger equations is solved using Duo with PECs of BeH + from Xu et al. 55 on a grid of 4001 points between 0.5 and 6.0 Å, evaluating 400 vibrational states (both bound and unbound) for X 1 S + , 393 for A 1 S + , 391 for B 1 P, and 381 for C 1 S + , imposing E max = 998 841 cm À1 . The X 1 S + state can hold 20 bound vibrational states with the rotational excitations ranging up to J max = 56. The photodissociation intensities were computed using Xu et al.'s transition dipole moment curves. For all temperatures, J max = 32 was used chosen as to maximize the agreement with the calculations of Yang et al. 56 Fig . 9 shows an example of how the photodissociation cross section of BeH + is recovered from the photoabsoption spectrum on top of the PECs involved. The bound-free transitions are given by dashed arrows, while the bound-bound transitions are depicted by the straight arrows. The inset in the figure shows the photoabsorption spectrum calculated at T = 1800 K. An estimate of the photodissociation cross section directly excluding the bound-bound contributions leads to the same trend as the previous method. Fig. 10 show the outcome of the two methods compared with the results from Yang et al. 56 in Panel A, the absolute difference between our two approaches is shown in Panel B. The greatest differences encountered at the boundaries, while the relative difference between the two methods is of the order of the 5 Â 10 À4 , making them equivalent for our purposes.
The photodissociation cross sections form two peaks with maxima at l = 172.31 nm and l = 96.01 nm. As the temperature increases, the height of the two peaks decreases alongside a flattening of the regions between the cross sections. In Fig. 9, our computed spectra (black) are compared with the values of Yang et al. 56 (red) for T = 1800 K, 4500 K, 10 000 K and 20 000 K. The photodissociation spectral contributions are obtained by subtracting the bound-bound components from the total photoabsorption spectrum, as shown in Fig. 9. The shapes and positions of the cross sections are the same in Yang et al. 56 and our model, but our model consistently gives lower cross sections (Fig. 11). This discrepancy is directly proportional to the temperature. An important contribution to the discrepancy between our results and those of Yang et al. 56 is the differences in the partition functions used. This difference is small for temperatures below 4500 K, where the two partition sums differ by about the 2% but is important for 10 000 K, where our partition function 10% bigger, and for 20 000 K, where our partition function is 30% bigger. Fig. 8 shows the difference of   55 and this work (DQ) as function of temperature, expressed in terms of percentage with respect to our (ExoCross) values. For temperatures below 7500 K, the difference between the two models is within the 2%, which increases at higher temperatures, with our model having higher values. The red dot points show the temperatures for which the photodissociation cross section are simulated. partition functions as function of the temperature; we believe ours to be more accurate (Fig. 11).

Conclusions
In this paper we develop a methodology which is suitable for the calculation of temperature-dependent photodissociation cross sections for diatomic molecules of arbitrary complexity. The methodology involves a solution of the coupled system of Schrödinger equations using bound vibrational basis functions and construction of pseudo-bound, temperature dependent cross-sections. For pure continuum cases, the photodissociation cross sections are generated using one of our smoothing approaches (Gaussian broadening or interpolation scheme). For hybrid, bound/free spectra, we use a stabilization method with averaged cross sections. We have tested the algorithm on a number of systems (HCl, HF, NaCl and BeH + ) containing both pure repulsive excited state potential curves and ones which support bound states. Agreement, within understandable limits, is found between published data and our calculations for all molecules examined. Our plan is to include theoretical temperature-dependent photodissociation cross sections of HCl, HF, NaCl and BeH + into the ExoMol data base, as part of the ExoMol project. 11 Results of these studies, after some improvements of the corresponding spectroscopic models, will be published elsewhere.
A major motivation for the developments presented here is the need to provide temperature-dependent photodissociation cross sections for polyatomic molecules such water and CO 2 . We note that the DVR3D program suite 57 has already been extended to provide rovibronic dipole transition intensities 58 and to use Lobatto shape functions 59 meaning that many of the developments for exploiting our proposed procedure are in place. However, for systems with more than two atoms, the identification of dissociative coordinates also becomes important. 59 We plan to extend our calculations of photodissociation cross sections and rates to triatomic and larger molecules. This will require us to develop a rigorous procedure appropriate for multichannel systems.

Data availability
Working examples of the input files used for generating the photodissociation cross sections and rates, together with the analysis tool are available as ESI † in the ESI.zip archive. A working version of Duo and ExoCross, which can be downloaded from the ExoMol github area, are required in order to run the input files. The file analysis.ipyn can be opened using Jupyter Notebook.
The directory duo-input contains the following input files that can be run using Duo, in order to create the initial state and transition files: 1 HCl-X-A.com, HF-X-A.com: input files for the A 1 P' X 1 S + band.
2 NaCl-Duo.com: input file for the A 1 P ' X 1 S + and the B 1 S + ' X 1 S + s bands.
The directory analysis contains the files required for calculating the cross sections and photodissociation rates: 1 hcl-2000.0.states and hcl-2000.0.trans: output from DUO calculation of HCl.
2 xsec-T2000.0.com: ExoCross input file needed for generating the raw absorption cross sections at T = 2000 K assuming the Gaussian profile of HWHM = 10 cm À 1 .

Conflicts of interest
There are no conflicts to declare.